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ABSTRACT 
The Joint Army/Navy Rotorcraft Analysis and Design (ANRAD) computer 


program was developed to aid in the analysis of helicopter rotor performance, stability 


and control, and rotor dynamics. JANRAD is an interactive, user friendly program, 


capable of accurately and quickly solving helicopter design problems at the 
preliminary design level. The program was written as a collection of MATLAB’ 
script and function files (M-files) using the 386-MATLAB’ version 3.5 programming 
language. The M-file janrad.m invokes the user interface routines and calls various 
analysis modules (M-files) which contain the appropriate analysis and output routines. 
Each of these modules use a common routine, trim.m, which employs blade element 
theory and a harmonic balance method for rotor trim. The program is limited to 


conditions of steady flight with no winds and is accurate at a hover and for forward 


airspeeds greater than or equal to 50 knots. 
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I. INTRODUCTION 
A. BACKGROUND | 

The development of rotorcraft design is a complicated iterative process involving 
the interaction of many variables. It is an especially difficult and time consuming 
task when the required calculations are performed by hand. There are several 
helicopter design software packages that perform the required calculations and 
automate the design analysis process. In searching for a computer program suitable 
for helicopter preliminary design, students found that current helicopter design codes 
were either too difficult to use, excessively detailed, o1 i:.:.uonally incomplete. As a 
result, the need for a custom computer code was identified. Through this thesis, the 
Joint Army/Navy Rotorcraft Analysis and Design (JANRAD) computer program was 
developed to meet the specific needs of the Naval Postgraduate School for preliminary 
helicopter design. 

JANRAD was written as an interactive, user friendly program, capable of 
accurately and quickly solving helicopter design problems at the preliminary design 
level. The program was intended Seiniaily for analysis of helicopter rotor 
performance, stability and control, and rotor dynamics. It was not intended to 
perform comprehensive helicopter design analyses. 

JANRAD was written as a collection of MATLAB’ script files (M-files). It was 
written specifically for version 3.5 of 386-MATLAB* [Ref. 1]. This was done for 


two reasons. First, MATLAB’ was ideally suited for the types of calculations and 


analyses involved in the preliminary helicopter design process. Second, engineering 
design students had typically used MATLAB’ for several projects and had a good 
working knowledge of the software and the programming language. The intent was to 
have a powerful and accurate computer aided design tool that one could master in a 
matter of minutes. In using JANRAD, students could occupy their time analyzing 
design issues rather than learning software. 

B. APPROACH TO THE PROBLEM 

The approach in developing JANRAD was to provide the design student with a 
rapid means of changing design parameters and obtaining accurate results over a wide 
range of helicopter configurations. To accomplish this, user input requirements were 
kept to a minimum, efficiency in parameter calculation routines was optimized, 
simplifying assumptions were made, and convergence criteria for iterative processes 
were carefully set to appropriate levels. 

To keep user input to a minimum, airframe and rotor geometric input 
parameters were minimized. The rotor blade was limited to a rectangular, uniform, 
rigid structure with linear twist. This allowed single value inputs for setup of the 
rotor blade parameters. The wing and tail section inputs were limited by using area 
and span to determine other geometric values. User input requirements were also 
reduced by requiring aircraft gross weight, equivalent flat plate area, auxiliary thrust, 
and aerodynamic surface coefficients of lift to be precalculated. 

To maximize the efficiency of the iterative processes, whenever possible, 


MATLAB’s’ built-in vector and matrix operations were used in place of repetitive 





"for" and “while loops. In addition, adjustments to the independent variables during 
the iterative processes were varied according to certain parameters such as advance 
ratio, density ratio, and magnitude of error. By varying the increment or decrement 
of the independent variables according to the appropriate parameters, the speed of 
convergence was greatly increased. 

In making assumptions and setting convergence criteria, the trade-offs between 
speed and accuracy were balanced by keeping in mind that JANRAD was developed 
for helicopter preliminary design. Accuracy was maintained for a wide range of 
helicopter configurations by using blade element theory to determine resultant forces 
and moments in the rotor system. Overall speed was maintained by making 
simplifying assumptions and selecting convergence criteria that resulted in small 
percentage errors that were acceptable for preliminary design. 

C. PROGRAM OUTLINE 

Currently, JANRAD consists of a group of 22 MATLAB’ script and function 
files. The main M-file, janrad.m, invokes the user interface functions and the various 
analysis modules. The routine allows the user to interactively input environment and 
design configuration parameters using the computer keyboard and screen. The user 
can either create or edit files containing parameters for specific designs and situations. 
The routine also provides the interface for selection of the various design analyses. A 
listing of janrad.m is found at Appendix A. Currently, there is a module available 
for rotor performance and one for stability and control. A rotor dynamics routine is 


under development. 


The heart of JANRAD is the M-file zrim.m. This M-file receives the 
environment and design configuration parameters from janrad.m and then balances the 
forces and moments in the rotor system for the appropriate airspeed. By trimming 
the forces in the rotor system, trim.m produces the parameters necessary for the 
various analysis modules. A listing of trim.m is found at Appendix B. 

The rotor performance module is called perf:m. This M-file accepts the data 
from trim.m and janrad.m, calculates the performance parameters, and generates the 
rotor performance output. Appendix C contains a listing of perf.m. The stability and 
control module, stab.m calculates and produces the stability and control output. Each 
M-file, with the exception of stab.m, is discussed in detail in this thesis. The details 
of the stability and control module are addressed in another thesis [Ref. 2]. 

Similarly, the details of the rotor dynamics module will be discussed in a future 


thesis. 














Il. PROGRAM DETAILS 

A. GENERAL 

JANRAD provides accurate data for standard or compound helicopters in 
hovering and forward flight. Forward flight analysis is only accurate for airspeeds of 
50 knots or greater. Conditions are limited to level flight with no winds. JANRAD 
provides data for one flight condition and one design configuration at a time. To 
obtain a range of data for an overall analysis, the user must make multiple runs of 
JANRAD. For each run, the appropriate parameters would be changed and the 
results could then be compared to those of previous runs. 
B. THEORY 

The underlying principle in helicopter rotor design is that the main rotor forces 
and moments must be adjusted so the rotor provides the lift and propulsive thrust 
required to meet the conditions of flight. The process of adjusting the rotor’s forces 
and moments is referred to as “trimming” the rotor. At a hover with no wind and 
assuming centered C.G. and no rolling moment due to the anti-torque device, there is 
no tilt required in the rotor disk. In this case, trimming the rotor involves adjusting 
collective pitch until the thrust generated by the rotor equals the weight of the aircraft 
plus the force created by the rotor downwash on the fuselage. In level, unaccelerated 
forward flight, assuming a centered C.G. and no rolling moments, the rotor disk must 
be tilted forward to create a vertical and a horizontal thrust component. The vertical 


thrust component must equal the weight of the aircraft (no downwash effects) minus 








lift created by any aerodynamic surfaces (assuming no lift generated by the fuselage). 


The horizontal component must equal the drag created by the fuselage and other 
aircraft components minus any auxiliary thrust. Trimming the rotor in forward flight 
involves the adjustment of cyclic pitch (collective pitch is the steady term of cyclic 
pitch) to maintain the required rotor disk tilt. 

For hovering and forward flight, JANRAD uses blade element theory to 
determine the rotor system forces and moments, and a method for harmonically 
balancing the forces and moments to adjust cyclic pitch [Ref. 3]. The only difference 
in the process for trimming the rotor system in hovering and forward flight is in how 
the rotor induced velocity is calculated. In forward flight, the rotor induced velocity 
is assumed uniform and is calculated using the inflow parameter A, [Ref. 4]. For 
hovering flight, induced velocity is determined by using an iterative process of 
equating the differential thrust calculated by momentum theory and blade element 
theory [Ref. 5, p. 38]. This process produces an elliptical non-uniform induced 
velocity distribution along the rotor blade. 

C. ASSUMPTIONS 

The following assumptions are made in order to reduce total computation time. 
The assumptions produce relatively small errors which are considered acceptable for 
helicopter preliminary design. 

1. Center of Gravity 

For the purposes of trimming the rotor system, the aircraft center of 


gravity is located at the center of the rotor hub. In doing so, the first-harmonic terms 





of thrust moment are eliminated. 
2. Lateral Forces and Moments 
The rotor is trimmed longitudinally and vertically. Lateral forces and 
moments due to the anti-torque device and aerodynamic surfaces are ignored. 
3. Rigid Rotor Blade 
For the purposes of the rotor trimming routine, the rotor blade is 
considered a rigid structure. This eliminates the need to consider blade modal 
responses and the subsequent affects on the rotor blade local angle of attack. 
4. Tip Loss 
Rotor blade tip loss is determined using a relationship involving coefficient 
of thrust. This results in a steady value of tip loss with respect to azimuth angle. 
The harmonic characteristics of tip loss are not considered. 
5. Aerodynamic Surfaces 
Lift and drag for the aircraft wing, horizontal tail, and vertical tail are 
calculated using two dimensional theory. Spanwise flow and tip losses are not 
accounted for. The aircraft fuselage lift characteristics are also ignored. 
6. Mach Number 
The airfoil lift and drag curves (coefficients of lift and drag versus angle of 
attack) for separate Mach numbers were reduced to one Mach independent curve. 
7. Wake Interactions 
Variations in angle of attack due to the interaction of rotor wake, 


aerodynamic surface trailing vortices, or fuselage interference were not accounted for. 





D. M-FILE DETAILS 
1. Janrad.m 

The main module for the JANRAD program is janrad.m. This M-file 
first queries the user to either edit or create a data file. If editing is selected, the 
program asks for the file name. If the file name is not found in the current directory, 
the program will tell the user to try again or instruct them on how to exit JANRAD. 
When the file name is found, a list of parameters is displayed, and the user is allowed 
to make changes as desired. If the user decides to create a new data file, the 
parameters are displayed one at a t*~2. In each instance, the user is prompted to 
input the desired value. Each prompt includes a display of the proper units. The 
parameters eligible for editing or creation are identical and are listed below: 


. pressure altitude 

. temperature 

airspeed 

aircraft gross weight 

. number of blades in the rotor system 

. rotor blade radius 

. rotor blade chord 

. effective hinge offset 

. rotor blade grip length 

0. rotor blade angle of twist 

11. rotor blade weight 

12. number of rotor blade elements 

13. rotor system rotational velocity 

14. number of rotor disk azimuth sectors 
15. average lift curve slope of rotor blade airfoil 
16. rotor blade airfoil 

17. collective pitch angle at 70% rotor blade radius 
18. equivalent aircraft flat plate area 

19. aircraft vertical projected area 

20. wing area 

21. wing span 

22. wing lift coefficient 


mOMIDWNAWNS 


23. wing profile drag coefficient 

24. wing efficiency factor 

25. horizontal tail area 

26. horizontal tail span 

27. horizontal tail lift coefficient 

28. horizontal tail profile drag coefficient 

29. vertical tail area 

30. vertical tail span 

31. vertical tail lift coefficient 

32. vertical tail profile drag coefficient 

33. auxiliary thrust 

Once the editing or creation process is complete, the program then queries 
the user for a file name. JANRAD saves the data to the current directory with the 
specified file name and a ".mat" extension. This is a binary file containing the file 
label and the names and values of the 33 parameters. 

Finally, the user is queried for a specific analysis. To date, the user can 
select rotor performance or aircraft stability and control. Janrad has provisions for a 
rotor system dynamics module. When the user selects rotor performance, the module 
perf.m is called. Similarly, for stability and control, stab.m is called. The user can 
also make selections to either change the current data (edit current parameters or 
create new ones) or exit the program. 

2. Trim.m 

The M-file trim.m trims the rotor system as appropriate for the 
environment and design parameters obtained as input in janrad.m. The rotor 
trimming process involves four steps: 
Calculate required parameters. 
. Guess rotor drag, cyclic pitch, and location of resultant thrust vector. 


. Trim the rotor system (calculate cyclic pitch). 
. Calculate rotor drag and location of resultant thrust vector. 


PWN 
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Initially, the program uses values for rotor drag, cyclic pitch, and the location of the 


resultant thrust vector based on established averages. The rotor system is then 
trimmed using those values and the other parameters previously calculated. Once the 
rotor is trimmed, rotor drag and the location of the resultant thrust vector are 
calculated. If the values for these two parameters fall within certain evaluation 
criteria (20% of prior rotor drag and 1.5% of prior location of resultant thrust 
vector), then the rotor is considered trimmed. If not, the rotor is retrimmed using 
new values for rotor drag and the location of the resultant thrust vector. This process 
is repeated until the two values fall within the evaluation criteria. 

As the first step, 7im.m calculates various parameters used throughout the 
routine. The parameters include ambient air density (9,), rotor drag or H-force 
(H,.o1), dynamic pressure (q), rotor disk area (Aji), rotor blade tip speed (V,,), 
fuselage drag (Dyuscuge), lift and drag on the aerodynamic surfaces (Dying, Dporiz wit» Dven 
wits Leings Lporiz wits Lven wi), the rotor tip path plane angle (aypp), and the advance ratio 
(u). The calculation of rotor drag at this point is a first guess based on the 


relationship: 


Ayocor = Va — (1) 


where: 


V.. = forward airspeed 


p, = ambient air density 


Pp, = sea level density 





The calculations for lift and drag on the wing and tail sections are based on a two 


dimensional analysis where: 


Ci (2) 

CoCo, * nmNeAR 
D=qC,S (3) 
L=qC,S (4) 


where: 

Cy, = drag coefficient 

C, = lift coefficient 

Cp, = profile drag coefficient 

e = wing efficiency factor 

AR = wing aspect ratio 

q = dynamic pressure 

S = wing area 
The values for the horiz:.tal and vertical tail efficiency factors are set at 0.8. The 
efficiency factor for the wing, as well as C,, and Cp, for the wing and other surfaces 
are direct input parameters obtained from janra:.m. The drag coefficients, aspect 
ratios, and areas of the aerodynamic surfaces, as well as q are calculated from inpt 


parameters obtained from janrad.m. The rotor tip path plan angle is calculated as: 
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D,+H 
@rpp= cans Pe Teesot 


Dy = Dyusetage + Dying + Doric wit + Deen at ~ Tux 
T,,, = auxiliary thrust 

Hao = H-force (rotor drag) 

GW = aircraft gross weight 

Ly = Lying + Lroriz wil 


Finally, the advance ratio is calculated as: 


_ ViCOS (@zpp) 
V, 


tip 

V.. = forward velocity 

Qypp = tip path plane angle 

Vip = rotor blade tip speed 

The required thrust and corresponding coefficient of thrust are also 


calculated. For hovering flight, thrust required is calculated as: 


rola (0 a Aus ps)| GW 
Aaisk 


Aven proj = Vertical projected area 


Ag = area of rotor disk 








GW = aircraft gross weight 


This equation assumes an effective vertical drag coefficient of 0.3 for all the aircraft 


components in the rotor wake [Ref. 5, p. 7]. In forward flight, the thrust required is 


calculated as: 
pense (8) 
COS (ppp) 
where: 
GW = aircraft gross weight 
L, = Liteg + Leora tit 
Qrpp = tip path plane angle 
The coefficiert of thrust is found by: 
Cy = et Ade =. (9) 


Agisk PaVeip 
where: 
T = thrust 
Aga = rotor disk area 
p, = ambient air density 


Vin = rotor blade tip speed 


Once C, is calculated, the rotor blade tip loss factor is determined as: 


=1-,| 2°r (10) 
Bei = 





where: 

C, = coefficient of thrust 

b = number of blades 

Once the required parameters are calculated, t7im.m sets up the rotor blade 
radial sections and rotor disk azimuth sectors. The radial sections are determined by 
evenly dividing the distance from the beginning of the aerodynamic portion of the 
rotor blade to the effective blade radius (R,, = B*R) by the number of blade elements 
input by the user in janrad.m. Azimuth is determined by dividing the 360 degree 
rotor disk into equally spaced areas determined by the number of azimuth sectors 
input by the user. 

The program then calculates rotor blade geometric angle due to twist (8) 
and rotor coning angle (8,). Blade twist is considered linear and determined by the 


relationship: 
B,=8,(0.7-) (11) 


Where: 

6, = blade twist 

r = distance from center of hub to center of blade element 

R = rotor blade radius 
Equation 11 takes this form because blade twist is converted to a positive value when 
input by the user in janrad.m. This allows the user to input either a positive or 


negative value of blade twist without affecting the accuracy of the program’s output. 
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It would be highly unlikely that blade twist would ever have positive value yet some 
may refer to it as positive when it is actually negative. The rotor coning angle is 


determined by: 


(12) 


where: 
6, = blade twist 
= distance from center of hub to center of blade element 


T = thrust 


b = number of blades 


R, = (Reg - ©)*1y 
R = rotor blade radius 
e = effective hinge offset 
rt; = location of resultant thrust vector 
Rig = effective radius (B*R) 
Woade = Weight of rotor blade 
Q = rotor rotational velocity 
Mise = Mass of rotor blade 
Equation 12 is based on a uniform rectangular blade, uses small angle approximations 


for 8,, and initially, assumes the location of the resultant thrust vector is at 0.7 r/R. 


15 


The next step is to determine the induced velocity distribution (v). Ata 


hover, v; is determined by equating differential thrust from momentum theory: 


dT=4p,nvirdr (13) 
and that of blade element theory: 
dT=p fs (a2)?a(6--2t) car (14) 
2 Qr 


which results in a quadratic equation of the form: 
2,/Q _{ Q2 a 
(4m) vj *(gbac|y, [Srpace| 0 (15) 


where: 
p, = ambient air density 


Q = rotor rotational velocity 


b = number of blades 
a = airfoil lift curve slope 
c = rotor blade chord 


6 = blade pitch angle (6, + 6.) 
8, = geometric angle due to twist 
|» = collective pitch angle 
r = distance from center of hub to center of blade element 
dr = blade element width 
The maximum root for v, is then taken and used in equation 14 to determine the 


differential thrust. These calculatiuns are repeated, incrementing 6 appropriately until 


16 





the sum of the dT’s equals the required thrust determined from equation 7. In 
forward flight, induced velocity is considered uniform and is calculated using the 


inflow parameter (A,;) where: 


1 





=C, 
Ap=p sin (arp,) + (16) 
Artp? 
The maximum root for \, is taken and v, is calculated as: 


where: 
». = advance ratio 
Qypp = tip path plane angle 
C, = coefficient of thrust 


Q 


rotor rotational velocity 

R = rotor blade radius 

V.. = forward airspeed 

The last step before starting the rotor trim routine is to make a first guess 
at cyclic pitch (8). Cyclic pitch is a first harmonic function with the cosine term (6,,) 
having a small positive value dependent on airspeed and the sine term (6,,), a negative 
value, also dependent on airspeed. The steady term (@,) or collective pitch term is an 


input parameter obtained from the user in janrad.m. The cyclic pitch is calculated as: 
6=8,+8,.cos () +6,,sin(y) (18) 
where: wy = azimuth angle 
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The rotor system is trimmed using blade element theory in a three step 
iterative process that consists of first adjusting collective pitch (@,), then trimming 
cyclic pitch (@), and finally retrimming the collective. Blade element theory is based 
on the determination of local angle of attack (a) for each rotor blade element. Once 
a is calculated, the coefficient of lift (C,) and coefficient of drag (C,) for each 
element are determined and finally, thrust and drag for each element are calculated. 


Blade element theory involves the use of the following equations: 


a =8+B.-6 (19) 
o-can[ (20) 
U, 
U,=v,cos (B,) +V.sin (@zpp) cos (B,) (21) 
+VCOS (@rpp) Sin (f,) cos (wy) 
U,=Qr+Vcos (a@,p,) sin (yp) (22) 


where: 
6 = cyclic pitch (eqn. 18) 
8, = rotor blade geometric angle 
¢ = inflow angle 
U, = velocity perpendicular to rotor blade 
U, = velocity parallel to rotor blade 
v; = induced velocity 


8, = coning angle 


18 





V.. = forward airspeed 
Q+pp = tip path plane angle 


yw = azimuth angle 





r = distance from center of hub to blade element 

The rotor trim procedure use two M-files, thrcalc.m and tmcalc.m. The 
listings for these M-files are found at Appendicies D and E respectively. The first M- 
file, thrcalc.m, uses blade element theory to calculate the differential thrust along the 
rotor blade at each azimuth angle. The differential thrust along the rotor blade is then 


summed to give the total thrust at each azimuth angle (T,). Differential thrust is 


calculated as: 


dT, = 4c dr (Up+U;) [C, cos (@) -C,sin(®) } (23) 


where: 
p, = ambient air density 
c = rotor blade chord 
dr = blade element width 


U, = velocity perpendicular to rotor blade 


U, = velocity parallel to rotor blade 


io) 
aS 
Ul 


coefficient of lift 
Cy, = coefficient of drag 
¢@ = inflow angle 


The small contribution of thrust provided by profile drag in the tip loss region of the 
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blade is accounted for by giving the profile drag coefficient (Cp,) a nominal value of 


.009. The equation for differential thrust in the tip loss region is as follows: 


aT ip" S PaC(R-Rese)(Up,,?+Uz,,2)(--0098in(b,4,)) (24) 


Peip 


where: 

p, = ambient air density 

c = rotor blade chord 

R = rotor blade radius 

R., = effective rotor radius 

U, up = velocity perpendicular to rotor blade 

Uniip = velocity parallel to rotor blade 

¢@ = inflow angle 
The second M-file, smcalc.m, also uses blade element theory to calculate the 
differential thrust moments along the blade at each azimuth angle. The differential 
thrust moments are then summed to give the total thrust moment at each azimuth 


angle (M,). The differential thrust moment is calculated as: 
dM, =p, cr dr (U;+ UZ) [C,cos () -C,sin() ] (25) 


where: 
p, = ambient air density 
c = rotor blade chord 


r = distance from center of hub to blade element 
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dr = blade element width 


Gq 
{ 


» = velocity perpendicular to rotor blade 


& 


= velocity parallel to rotor blade 
C, = coefficient of lift 
Cp = coefficient of drag 
¢@ = inflow angle 
The small thrust moment provided by profile drag in the tip loss region of the blade is 


accounted for in a similar manner as the thrust in thrcalc.m. The equation is: 


rt (R-Ry ge) (26) 
2 


OM. p= 3 pac(r (R-Reee) (Up? + Ur...) (-.009 sin (dip) ) 


p, = ambient air density 

c = rotor blade chord 

R = rotor blade radius 

Rig = effective rotor radius 

U, ip = velocity perpendicular to rotor blade 

Uiip = velocity parallel to rotor blade 

¢@ = inflow angle 

The two M-files, thrcalc.m and tmcalc.m, call an M-file which contains the 
lift coefficient and the drag coefficient data for the appropriate rotor blade airfoil. 
The airfoil data are arranged as separate polynomials for various ranges of angle of 


attack (a). There are currently two airfoil data function files, hhO2clcd.m and 
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vrl2clcd.m. The former contains data for the McDonnell Douglas Helicopter 
Corporation (MDHC) HH-02 airfoil used on the AH-64A Apache attack helicopter. 
The latter contains data for the Boeing VR-12 airfoil planned for use on the RAH-66 
Comanche. There is no Mach number dependency incorporated in either airfoil data 
function files. The listings for these M-files are found in Appendicies F and G 
respectively. 

In both the initial trim and retrim portions of the collective trim procedure, 
the collective is adjusted by increasing or decreasing @, a suitable amount to match the 
required thrust calculated from either equation 7 or 8. The cyclic terms (6,, and 6,,) 
are not adjusted during the collective trim procedure. Collective pitch is adjusted 
until the resulting thrust calculated by thrcaic.m is within 2% of required thrust for 
the initial trim procedure and 1% of required thrust for the retrim procedure. 

The cyclic trim routine begins by determining the initial thrust moment 
(Mi) based on the cyclic pitch determined by the collective trim procedure. This 
cyclic pitch value is considered the initial pitch (6). The thrust moment is 
calculated using the M-file yncalc.m. An initial partial derivative of cyclic pitch with 
respect to thrust moment is determined by collectively increasing cyclic pitch a very 


small amount and determining the corresponding thrust moment. Giving the new 


cyclic pitch and corresponding new thrust moment the labels 0,44, and Maew 


respectively, the derivative is determined as: 








(new) ~ 


6 8 ine) (27) 


RB. 
OM Minew ~ Meine) 


Next, the first harmonic terms are removed from the original thrust moment (M,,,). 
Using a numerical Fourier analysis, the first harmonic terms of the thrust moment 


are: 


n 
M, = 2 Minow Sin (W) (29) 
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where: 
n = number of azimuth sectors 
M, = total thrust moment at each azimuth angle 
y = azimuth angle 


The new thrust moment (M,,,.)) is then calculated as follows: 
Minew = Minty 7 i- COS (W) -M,,sin (py) (30) 
The change in thrust moment is then calculated as; 
AM=M new ~ Mine) (31) 


and the new cyclic pitch (6,,..)) is then determined by the relationship: 


008 
8 (new) = G (ine) > vad (32) 


A first harmonic function is then enforced on theta by using equation 18 where: 
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n 
6, .= > 8 new COS (HW) 


n 
aD 8 new 8in (W) (34) 


Using the new cyclic pitch from equation 32, a new thrust moment is calculated, and 
subsequently a new change in moment is determined. The previous values for cyclic 
pitch and thrust moment become the initial values, and the process is repeated 
(skipping the small collective increase used for the initial determination of the 
derivative) until there is only a small first harmonic (2% of the total thrust moment) 
left in the thrust moment. 

At anytime during the collective or cyclic trimming routines, if the values 
of the respective parameters (difference in thrust or difference in thrust moment) 
begin to diverge, program execution will halt. Divergence of these parameters 
indicate that the blade is stalling over too large of an azimuth region. At this point, 
the user is told to use a slower airspeed or redesign the system. 

Once the cyclic is trimmed and the collective retrimmed, trim.m, using 
dmcalc.m, then calculates the differential drag values for each azimuth angle (dD,) 
and the rotor drag moments for each azimuth angle (DM,). A listing for dmcalc.m is 
found at Appendix H. The drag moments are determined by summing the product of 
the differential drag values and radial distance from the center of the hub to the center 
of each differential blade element. The differential drag values are calculated as 


follows: 





dD, = +p, Cdr(Up+U2) (C,sin (o) + Cycos (6) ] (35) 


where: 

p, = ambient air density 

c = rotor blade chord 

dr = blade element width 

» = velocity perpendicular to rotor blade 

U, = velocity parallel to rotor blade 

C, = coefficient of lift 

Cy = coefficient of drag 

¢@ = inflow angle 
The contribution of profile drag in the tip loss region is accounted for in the drag and 
drag moment equations in a similar manner as in thrcalc.m. The equation for the 


differential drag in the tip loss region is: 
1 
OD. ip= SPac (R-Re ge) (Up.ie + U0.) (.009cos (9,,;,) ) (36) 


where: 
p, = ambient air density 
c = rotor blade chord 
R = rotor blade radius 
R.y = effective rotor blade radius 


U, up = velocity perpendicular to rotor blade 
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U.u, = velocity parallel to rotor blade 


¢i» = inflow angle 


The equation for the drag moment in the tip loss region is: 


aDM, 


tip 


(R-Roee) 
2,4 R- {eRe (37) 


Once the differential drag values are determined, the H-force (rotor drag) 
and the location of the resultant thrust vector are then calculated. The H-force is a 


first harmonic function and is calculated as follows: 


bcos (a,,,) (< ; R 
Hyoror = Pees eee) (S H,,- sin (B,) > Hoe (38) 
Zo to 


where: 
= number of blades 
Qrpp = tip path plan angle 
r, = Start of aerodynamic portion of blade 
R = radius of rotor blade 
8B, = coning angle 


The cosine and sine terms of the rotor H-force are determined as: 


n 

Me" Se aT, cos (wy) (39) 
is 

y= =) dp, sin() (40) 


i*1 


If the difference between the initial rotor H-force and the computed value is more 





than 20% of the initial one, then the location of the resultant thrust vector is 
computed and the rotor trim process is repeated. If the difference is less than 20%, 
then the location of the resultant thrust vector is computed and the rotor trim process 


is repeated only if the new location of the resultant thrust vector does not fall within 


its limits. 
The location of the resultant thrust vector (r;) is calculated as follows: 
Lee mat (41) 
Ty 
where: 


M, = mean of the sum of thrust moments 

T, = mean of the sum of thrust 
If the difference between the initial r; and the computed one is more than 1.5% of the 
original, then the rotor trim process is repeated. If the difference is less than 1.5%, 
then the rotor is considered trimmed. 

3. Perf.m 

This M-file calls the rotor trim routine, calculates various rotor 

performance parameters, and then displays and saves the rotor performance data. 


After calling trim.m, the routine calculates rotor power required and the resulting 


rotor torque using the relationships: 


EDM 0 


nD (42) 


Protor (h.D- ) “S50 
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- 2M 


rotor = 


DM, = total drag moment at each azimuth angle 

n = number azimuth sectors 

b = number of blades 

Q = rotor rotational velocity 

The routine then calculates rotor solidity (0), coefficient of thrust over 


solidity (C,/0), coefficient of torque over solidity (Cg/o), coefficient of H-force over 


solidity (C,,/a), disk loading (D.L.), figure of merit (F.M.), and advancing blade tip 


Mach number (M,,). Solidity is calculated as: 


b = number of blades 
c = blade chord 
R = blade radius 


C, was previously calculated by trim.m and Cx and Cy, are calculated as: 











= Orotor (45) 


Co= 
Aaisx Pa Veip R 


H. 
Cy= rotor (46) 


2 
Aagisk Pa Veip 


where: 


Qurosor 


Agi = area of rotor disk 


rotor torque 


p, = ambient air density 


Vip = blade tip speed 


Huo, = rotor H-force 
C;, Cg, and C,, are then divided by o to get the C,/o relationships. Disk loading and 


figure of merit are determined by: 


Peli, (47) 
714 
D.L. 
Tyotor 20, (48) 
F.M.= 
550 (Protor) 


where: 
Toor = rotor thrust 
R = rotor radius 
p, = ambient air density 
P cor = rotor power (h.p.) 


Disk loading and figure of merit are set to zero for hovering flight. The advancing 
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tip Mach number is calculated using the relationship: 


2 Ve ip COS (Grp) +VL (49) 


Mee A9, 05 ,/T,+460 
where: 

Vip = rotor blade tip speed 

Qrpp = tip path plane angle 

V.. = airspeed 

T, = ambient temperature (°F) 

Once the rotor performance parameters are calculated, perf.m then asks the 
user if they want the rotor performance parameters displayed on the screen. If the 
answer is yes, the data is presented as two screens of input data and two of calculated 
data with a user controlled pause between screens. 

Finally, perf.m saves the rotor performance data to disk as two separate 
files. The first file contains the single valued input and calculated parameters and is 
given the file name the user defined in janrad.m with a " prf" extension. This file is a 
MS-DOS" text file. The file is formed using the MATLAB* DIARY function which 
records a series of MATLAB’ FPRINTF commands and then stores them to disk as a 
file named diary. A call to the operating system is then made renaming diary to the 
appropriate ".prf" file name. The second file contains the vector si matrix valued 
data and is given the user defined file name concatenated with an "_p" and a ".mat" 
extension. This file is a binary file which the user can access using the MATLAB” 


LOAD command. The file is created using the MATLAB® SAVE command. 
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4. Stab.m 
The M-file, stab.m, performs the helicopter stability and control 
calculations, formats the stability and control output, and saves that output to disk. 
This module uses the following M-files: 1) cbodygrp.m; 2) cmdbwplc.m; 3) 
cmdbwplh.m; 4) cmrgrp.m; 5) cruise.m; 6) ctrgrp.m; 7) dctplots.m; 8) dctplott.m; 9) 
deriv2.m; 10) hmrgrp.m; 11) hover.m; 12) hirgrp.m; and 13) stabout.m. For details 


on stab.m and the various routines, See reference 2. 
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Iii. USER INSTRUCTIONS 
A. GENERAL 
Program installation, guidance for input requirements, and rotor performance 
output are discussed in the subsequent paragraphs. These instructions assume that the 
reader has a working knowledge of the MS-DOS* operating system and of 386- 
MATLAB’ version 3.5. 
B. INSTALLATION 
To install JANRAD: 
1. run MATLAB’ 
2. ensure you are in your working directory 
3. from the MATLAB* command prompt, type »!drive: (A: or B:) then 
type »INSTALL. The M-file, install.m, will create a subdirectory named JANRAD. 
It will then copy the JANRAD files to that directory. Once the files are copied, the 
current directory is changed to JANRAD and the program is ready to run. 
C. PROGRAM EXECUTION 


1. Input 


To run JANRAD, at the MATLAB’ command prompt, type »JANRAD. 


At the initial screen, enter a 1 if you want to edit a previously stored data file, 
otherwise enter a 2. If you chose to edit a file, you will be prompted to enter the file 
name without the extension. If JANRAD cannot find the file you entered in the 


current directory, it will allow you to try again. Check your spelling or insure you 
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are in the JANRAD directory. Once the file is loaded, the edit menu, as seen in 


Figure 1, will be displayed. 


*** EDIT MENU *** 


. pressure altitude . temperature 

. airspeed . gross weight 

. number of blades . blade radius 

. blade chord . hinge offset 

. blade grip length . blade twist 

. blade weight . # blade elements 
rotational velocity . # azimuth sectors 
lift curve slope . airfoil 

. collective pitch . flat plate area 

. vert projected area . wing area 


. wing span . wing CL 

. wing CDo . wing efficiency factor 
. horizontal tail area . horizontal tail span 

. horizontal tail CL . horizontal tail CDo 

. vertical tail area . vertical tail span 

. vertical tail CL . vertical tail CDo 

. auxiliary thrust 


- NO CHANGES 


Input the parameter to change: 





Figure 1. Edit Menu 


Select the number of the parameter you wish to change. The current value for that 
parameter will be displayed followed by a prompt to enter the new value. If you 
choose not to change the value, simply press <Enter> and the previous value will be 
maintained. Once a value is entered, or <Enter> is pressed, the editing menu will 
again be displayed. Choose the number of the next parameter to change, or a 0 to 
exit the editing menu. If you chose to create a new data file, prompts for individual 
parameters will be displayed. Enter the value for each parameter when prompted. 

Enter values for the specified parameters regardless of editing or creating a 
data file in the following manner: 


1. PA - Enter the pressure altitude in feet. 
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2. temperature - Enter the temperature in degrees fahrenheit. 

3. airspeed - Enter forward airspeed in knots. JANRAD automatically converts 
airspeed to feet/second. 

4. GW - Enter aircraft gross weight in pounds. 

5. number of blades - Enter the number of blades used in the rotor system. 

6. radius - Enter the rotor blade radius in feet. Measure the radius from the 
center of the rotor hub to ti.e tip of the blade. 

7. chord - Enter the average rotor blade chord in feet. 

8. hinge offset - Enter the effective hinge offset in feet. 

9. blade grip - Enter the distance from the center of the rotor hub to the 
beginning of the aerodynamic portion of the rotor blade in feet. 


10. blade twist - Enter rotor blade twist in degrees. You can enter blade twist 


as either a positive or negative value. JANRAD uses the absolute value of blade 


twist. 
11. blade weight - Enter the weight of a single rotor blade in pounds. 
12. rotor blade elements - Enter the number of rotor blade radial elements in 
. JANRAD adds one element to account for tip loss. 
13. rotational velocity - Enter the rotor rotational velocity in radians/second. 
14. azimuth sectors - Enter the number of rotor disk azimuth sectors. 
15. lift curve slope - Enter the average slope of the linear portion of the rotor 
blade airfoil’s lift curve, where the lift curve is C, versus alpha. 


16. rotor blade airfoil - Enter a 1 if you are using an HH-02 airfoil. Enter a 2 





if you are using a VR-12 airfoil. The HH-02 airfoil data is very close to the NACA 
0012 airfoil data. Therefore, enter a 1 for a generic airfoil. 

17. collective pitch - Enter the estimated collective pitch at 0.7 r/R in degrees. 
JANRAD automatically converts the value to radians. 

18. flat plate area - Enter the aircraft equivalent flat plate area in square feet. 

19. vertical projected area - Enter the vertical projected area in square feet. 

20. wing area - Enter the wing area in square feet. Enter a 0 if the aircraft has 
no wing. 

21. wing span - Enter the wing span in feet. If you included the segment of the 
fuselage encompassed by the wing in the wing area value, include it in the span value 
also. 

. wing C, - Enter the expected lift coefficient for the wing. 

. wing Cp, - Enter the wing profile drag coefficient. 

. wing efficiency factor - Enter the wing efficiency factor. 

. horizontal tail area - Enter the horizontal tail area in square feet. 


. horizontal tail span - Enter the horizontal tail span in feet. 


. horizontal tail C, - Enter the expected lift coefficient for the horizontal tail. 


. horizontal tail Cp, - Enter the horizontal tail profile drag coefficient. 

. vertical tail area - Enter the area of the vertical tail in square feet. 

. vertical tail span - Enter the span of the vertical tail in feet. 

. vertical tail C, - Enter the expected lift coefficient for the vertical tail. 


. vertical tail Cp, - Enter the vertical tail profile drag coefficient. 
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33. auxiliary thrust - Enter the expected auxiliary thrust in pounds. 


Once the edited or new values are entered, you will be prompted to save 


the data under a designated file name. The file name must use a combination of six 


letters, numbers, or the underscore. If more than six characters are entered, you will 
be prompted to try again. JANRAD will save the file to the current directory 
(JANRAD) as the file name with a ".mat" extension. If you were editing a data file, 
press <Enter> if you want to save the file with the original file name. 

The execution menu is displayed next. Enter a 1 if you want to perform a 
rotor performance analysis. Enter a 2 if you choose to perform a stability and control 
analysis. Enter a 3 to execute the rotor dynamics analysis. The rotor dynamics 
module is currently under development. If you enter a 3, this fact will be displayed 
and you will be asked to enter another selection. If you decide to edit your data 
further or you want to create a new file, enter a 4. To quit JANRAD, enter a 5. 

From any of the selection menus, if you enter a character other than what 
is requested, JANRAD will request that you enter an appropriate response. If you 
enter a letter or special ASCII character, the MATLAB’ prompt “undefined variable 
or function" will be displayed along with a repeat of the initial prompt. If you enter a 
number other than requested, the initial prompt will be repeated. 

No further input is required if you select rotor performance. See 
Reference 2 for details on the input required if you select stability and control. 


Details for the rotor dynamics input requirements should be covered in a future thesis. 








2. Output 


If rotor performance is chosen from the execution menu, you will be 


asked if you want the rotor performance results displayed on the screen. If you 


respond in the affirmative, four screens will be displayed. The screens will look 


similar to those shown in Figures 2 through 5. The data in these figures is from the 


example helicopter (Prouty’. helicopter) found in Reference 5, page 669. The first 


two screens contain the input data and the second two contain the output data. 


*«* INPUT DATA *** 


prouty 


Forward velocity 
Temperature 
Pressure altitude 
Gross weight 
Number of blades 
Rotor radius 

Blade chord 

Blade twist 

Blade lift curve slope 
Blade weight 
Rotational velocity 
Blade grip length 
Hinge offset 


press any key to continue... 


Figure 2. First Output Screen (Input Data) 


= 
= 
= 
= 
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110 
59 





kts 
degs F 
ft 
lbs 


ft 
ft 
degs 


lbs 
rads/sec 
ft 

ft 


*** INPUT DATA CONTINUED *** 
prouty 


f£t*2 
ft*2 
ft*2 
ft 


Equivalent flat plate area 
Vertical projected area 
Wing area 

Wing span 

Wing CL 

Wing CDo 

Wing efficiency factor 
Horizontal tail area 
Horizontal tail span 
Horizontal tail CL 
Horizontal tail CDo 
Vertical tail area 
Vertical tail span 
Vertical tail CL 
Vertical tail CDo 
Auxiliary thrust 


ft*2 
ft 


press any key to continue... 





Figure 3. Second Output Screen (Input Data Continued) 


*** CALCULATED DATA *** 
prouty 


1232 
235 
0 

0 
+222 
3 
813 
114 


Fuselage drag 

Rotor drag 

Wing lift 

Wing drag 

Horizontal tail lift 
Horizontal tail drag 
Vertical tail side force 
Vertical tail drag 


-48 
-69 
-63 
-98 
.07 
32 


Tip path angle 

Rotor coning angle 

Location of mean thrust (r/R) 
collective pitch at .7 r/R 
lst lat cyclic term-Al (deg) 
lst long cyclic term-Bl1 (deg) 


press any key to continue... 





Figure 4. Third Output Screen (Calculated Data) 
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*** CALCULATED DATA CONTINUED *** 
prouty 


0.085 
0.00 lbs/ft*2 
0.00 
0.084 
.0034 
.0010 
0.746 
0.285 
20284 lbs 
959 h.p. 
24347 ft-lbs 


solidity 

Disk loading 
Figure of Merit 
CT/sigma 


z 
= 
= 
= 
CQ/sigma = 0 
CH/sigma = 0 
Tip mach of the adv. blade = 
Advance ratio = 
Rotor thrust required (TPP) = 
Rotor power required = 
Rotor torque = 





Figure 5. Fourth Output Screen (Calculated Data Continued) 
Each of these screens is displayed only one time. Once you press any key to go to 
the next screen, you cannot return to the previous one. A copy of the four screens is 
saved to the current (VANRAD) directory with the file name you entered at the “save 
instructions" screen and concatenated with a “.prf" extension. 

The vector and matrix data generated by the rotor performance module of 
JANRAD is saved to the current (JANRAD) directory with the file name you entered 
at the “save instructions” screen and concatenated with a "_p" and a ".mat" extension. 
The following vector and matrix parameters are saved: 

1. r- a vector containing the radii (in feet) to the center of each blade element. 
This is a row vector containing n elements, where n = # of blade elements + 1. 

2. yw =a vector containing the azimuth angles (in degrees) around the rotor 
disk. This is a column vector containing n elements, where n = # of azimuth 
sectors. 

3. v,; = a vector containing the induced velocity distribution (in ft/sec) along the 


rotor blade. This is a row vector containing n elements, where n = # of blade 
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elements + 1. 

4. 6 = vector containing cyclic pitch (in degrees) with respect to azimuth 
angle at 0.7 r/R. This is a column vector containing n elements, where n = # of 
azimuth sectors. 

5. 8, = a vector containing blade pitch (in degrees) along the rotor blade. This 
is a row vector containing n elements, where n = # of blade elements + 1. 

6. a@ = a matrix containing the angle of attack (in degrees) distribution. This is 
an m Xn matrix, where m = # of azimuth sectors and n = # of blade elements +1. 

7. Ty, = a vector containing the total thrust (in pounds) at each azimuth sector. 
This is a column vector containing n elements, where n = # of azimuth sectors. 

8. M, = a vector containing the total thrust moment (in ft-lbs) at each azimuth 
sector. This is a column vector containing n elements, where n = # of azimuth 
sectors. 

9. DM, = a vector containing the total drag moment (in ft-lbs) at each 
azimuth sector. This is a column vector containing n elements, where n = # of 
azimuth sectors. 

10. dT = a matrix containing the differential thrust (in pounds) distribution. 
This is an m Xn matrix, where m = # of azimuth sectors and n = # of blade 
elements +1. 

11. dM = a matrix containing the differential thrust moment (in ft-lbs) 
distribution. This is an m Xn matrix, where m = # of azimuth sectors and n = # of 


blade elements +1. 





12. dD = a matrix containing the differential drag (in pounds) distribution. 


This is an m Xn matrix, where m = # of azimuth sectors and n = # of blade 





elements +1. 

A final rotor performance module screen is displayed instructing you on 
how to access your rotor performance data. Use the MS-DOS° TYPE command or a 
DOS editor to view or print the single variable data. The vector and matrix data can 
be loaded into memory using the MATLAB’ LOAD command. Some example plots, 
using vector and matrix data from the Prouty helicopter, are seen below. Figure 6 is 


a plot of hover induced velocity (ft/sec) versus fraction of blade radius (r/R). 





Figure 6. Hover Induced Velocity vs. Blade Radius 





A plot of total thrust moment (ft-lbs) versus azimuth angle (degrees) for the Prouty 


helicopter at 125 knots is shown in Figure 7. 





Figure 7. Total Thrust Moment vs. Azimuth Angle 


A three dimensional plot showing the differential thrust (Ibs) distribution over blade 
radius versus azimuth angle is at Figure 8. Figure 9 shows the angle of attack 
(degrees) distribution. The data for both plots is from the Prouty helicopter at 125 


knots. 
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Figure 8. Differential Thrust Distri 


Figure 9. Angle of Attack Distribution 
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The plots in Figures 6 through 9 were generated from the data contained in a data file 
named prouty_p.mat. The data was plotted using MATLAB” 
for WINDOWS’, version 4.0. 

The preceding paragraphs discussed the details for the rotor performance 
output, the details for stability and control output are contained in Reference 2. The 
details for the blade dynamics module will be found in a future thesis. 

D. PROGRAM UPDATES 

JANRAD is written entirely in the 386-MATLAB’ version 3.5 programming 
language. This version of JANRAD is designated version 1.0. The text in this 
document corresponds to the software loaded on a master 3.5" diskette labeled 
JANRAD version 1.0, September 1993. Subsequent changes to the software should be 
given a new numerical designation and loaded on a separate diskette. Each new 
version of MATLAB* utilizes a slightly different command set. As a result, slight 
modifications may be required to run JANRAD version 1.0 on later versions of 
MATLAB’. Additionally, JANRAD version 1.0 was developed for use on computers 
using the DOS type operating system. JANRAD makes some calls to the operating 
system using the MATLAB* "!" function. For this reason and because actual 
MATLAB’ commands differ slightly, modifications would be required to run 
JANRAD version 1.0 on an operating system other than (i.e., UNIX® ). The ".mat" 
data files created by JANRAD version 1.0 can be transferred to other operating 
systems using a file transfer program. Those files could then be loaded by the 


particular operating system’s version of MATLAB’ and used for plotting. 


IV. RECOMMENDATIONS 

Improvements to this program code can be made in several areas. First, the 
blade dynamics module should be developed as a follow-on thesis or series of AE 
3402 (Helicopter Dynamics) class projects. As another thesis project, the program 
code should be rewritten to work with the UNIX® version of MATLAB’. The UNIX” 
version runs faster and has enhanced graphic capabilities. The program code should 
also be updated to include: 

1. The capability to offset the center of gravity from the center of the rotor hub 
during the rotor trimming process. 

2. A method of incorporating Mach number dependency in the rotor blade 
airfoil data. This may require the use of "look-up" tables versus the current use of 
parametric equations. 

3. A more accurate method for determining rotor induced velocity in forward 
flight. This would be especially beneficial for airspeeds less than 100 knots. 

4. Changes to the code to account for the effects of rotor downwash on the 
wing angle of attack. 

5. Changes to the code to incorporate a non-rectangular, non-uniform, flexible 
rotor blade. 

These updates to the code would expand the design students analysis capabilities as 
well as improve the overall accuracy of the program. This would be especially true 


when utilizing the enhanced capabilities of the UNIX® version of MATLAB’. 
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It is also recommended that this version of JANRAD (version 1.0) be kept on a 
master diskette and that future updates be stored on a separate disk. This will ensure 
that a working copy of the program is always available. Thesis level updates to 
JANRAD should, at a minimum, refer to this document. Updates that are not large 
enough to justify thesis work should be documented as attachments to this document. 
It is recommended that all update documentation be maintained (chronologically) in a 


binder along with this document. 
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APPENDIX A: 


JANRAD .M 


Joint Army Navy Rotorcraft Analysis and Design 
(JANRAD) 


Version 1.0 
September 1993 


MAJ Bob Nicholson 
MAJ Walter Wirth 


This program is an interactive preliminary design tool 
developed to aid the design student in determination of 
initial rotorcraft configurations and in the calculation 

of performance, stability and control, and other parameters. 
The program will work for conventional or compound rotorcraft. 
It will provide accurate data for airspeeds less than 10 


knots and greater than or equal to 50 knots. 


Variable List for janrad.m, 


trim.m, thrcealce.m, tmcalc.m, dmcalc.m, 


hhO2cled.m, vrl2clcd.m, and perf.m 


a lift curve slope of rotor system airfoil 
Adisk area of rotor disk 
Afh fuselage equivalent flat plate drag area 
Afv vertical projected area (fuselage area under disk) 
airfoil rotor system airfoil type (HH02/VR12) 
alpha angle of attack, rotor blade radial segment 
alphaT rotor tip path plane angle 
b number of rotor blades 
B tip loss parameter 
betao rotor coning angle 
betat geometric angle, rotor blade radial segment 
bhoriz span, horizontal tail 
bvert span, vertical tail 
bwing span, wing 
cblade chord, rotor blade 
cD drag coefficient, rotor blade radial segment 
CDohoriz profile drag coefficient, horizontal tail 
CDovert profile drag coefficient, vertical tail 
CDowing profile drag coefficient, wing 
CDhoriz drag coefficient, horizontal tail 
CDvert drag coefficient, vertical tail 
CDwing drag coefficient, wing 
CH rotor H-force coefficient 
CH_sig CH/solidity 
CL lift coefficient, rotor blade radial segment 
CLhoriz lift coefficient, horizontal tail 
CLvert lift coefficient, vertical tail 
CLwing lift coefficient, wing 
cQ rotor torque coefficient 
CQ_sig CQ/solidity 

rotor thrust coefficient 
CT_sig CT/solidity 
aD differential drag, rotor blade radial segment 
dadD differential drag, rotor blade tip 
adDM differential drag moment, rotor blade tip 
daM differential thrust moment, rotor blade tip 
adT differential thrust, rotor blade tip 
de1lM change in total thrust moment 
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Dftotal 
Dfuse 
DL 

aM 


DMpsi 
dr 


Drotor 


Hrotor 
lamdaT 
Lftotal 
Lhoriz 
Lvert 
Lwing 
Mic 
Mis 
Machtip 
mblade 
Mpsi 
mu 

naz 

nbe 
omega 
PA 

phi 
phitip 
Protor 
psi 


Qrotor 
x 

R 

Rbar 
RbarT 
Reff 
rho 

rT 
solidity 
Shoriz 
Svert 
Swing 
T 

Taux 
temp 
theta 
thetalc 


resultant of fuselage drag and aux thrust 

total drag generated by non-rotor bodies 

disk loading 

differential thrust moment, rotor blade radial seg 
total blade drag moment at specific azimuth angle 
rotor blade radial segment width 

rotor system drag 

differential thrust, rotor blade radial segment 
drag, horizontal tail 

change in cyclic pitch with change in thrust moment 
drag, vertical tail 

drag, wing 

effective hinge offset 

wing efficiency factor 

name of input file 

figure of merit 

length of inner non-aerodynamic portion of blade 
aircraft gross weight 

rotor H-force 

forward flight induced velocity parameter 

total lift generated by non-rotor bodies 

lift, horizontal tail 

lift, vertical tail 

lift, wing 

first harmonic (cosine) thrust moment coefficient 
first harmonic (sine) thrust moment coefficient 
Mach number at rotor blade tip 

mass of rotor blade 

total blade thrust moment at specific azimuth angle 
advance ratio 

number of azimuth sectors 

number of blade elements 

rotor rotational velocity 

pressure altitude 

inflow angle, rotor blade radial segment 

inflow angle, rotor blade tip 

power required by rotor 

azimuth angle 

dynamic pressure 

rotor torque 

radius, rotor blade radial segment 

rotor blade radius 

Reff-e 

rT*Rbar 

effective rotor blade radius (tip loss) 

ambient air density 

location of resultant thrust vector 

solidity 

area, horizontal tail 

area, vertical tail 

area, wing 

rotor thrust 

auxiliary thrust 

ambient air temperature 

cyclic pitch 

first harmonic (cosine) of cyclic pitch 

first harmonic (sine) of cyclic pitch 
collective pitch at .7 r/R 

total blade thrust at specific azimuth angle 
geometric rotor blade twist 

vertical component of velocity 

vertical component of velocity at tip 


49 








* Ut horizontal component of velocity 

*  Uttip horizontal component of velocity at tip 

% vi induced velocity 

* Vinf forward airspeed 

*  Vtip tip speed 

% wblade weight of rotor blade 

cle 

disp(’ ‘) 

disp(’' ‘) 

disp (’ TRA AASZESAALELESZEEERESS EE EE SETTER EEE PES CCST ECS ECS CSC CASES SLES ALS ) 
disp (’ * “t) 
disp (’ * *** Joint Army/Navy Rotorcraft Analysis and Design *** *) 
disp (’ . (JANRAD) **) 
disp (’ * mie) 
disp (' * Version 1.0 wr) 
disp (’ * September 1993 #e) 
disp (’ * *?) 
disp (’ PESSRERELALALE SELES EERE PERE LES EPPS TET TC CSS CPCS CSCC SCC SC CEL PCP eee es ) 
Gisp(’ ‘) 

checkl=1; 


while checkl > 0 
* clearing all variables in the MATLAB environment 


clear 
checkl=1; 


disp(’ ’) 

disp (’ Do you want to edit an existing file or create a new one?’) 
check=1; 

while check > 0 


answer=input (’ 1. edit existing file 2. create new file >> '); 


% «** If editing an existing file: get file name, display edit 
% menu, allow changes to selected variables, and save under 
% desired file name. Loads to and saves from current 

% directory as a .mat file. *** 


if answer==1, 


cle 
disp(’ ') 
disp(’ ’) 
disp (’ *** LOAD INSTRUCTIONS *** ‘) 
disp(’ ‘) 
disp (’ A. Input the name of the file to edit.’) 
disp (’ B. The file was saved in your previous session’) 
disp (’ with a ".mat" extension.’) 
disp (‘ C. Do not include the extension or quotations.’) 
disp(’ °) 
disp (’ ex: dsgnl’) 

flag=0; 

while flag < 1 

filenamel=input (' name of input file: ‘','s’); 
eval (['flag=exist (''’,filenamel,’.mat’’);’)) 
if flag < 1, 
disp(’ ‘) 
disp (‘ The file does not exist, try again or < Ctrl-C >’) 
disp (’ to exit program.’ ) 
end 
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end 
eval({‘load ’,filename2)]) 


disp (' 
disp (’ 
Gisp (’ 
disp (’ 
disp (’ 
disp (' 
disp (’ 
disp (’ 
disp (’ 
disp (’ 
disp (’ 
disp (' 
disp (' 
disp (’ 
disp (' 
disp (’ 
disp (’ 
disp (’ 
disp (’ 
disp (’ 
disp (’ 
disp (’ 


while check > 0 


ae 
«** EDIT MENU 
") 
1. pressure altitude 
3. airspeed 
5. number of blades 
7. blade chord 
9. blade grip length 
11. blade weight 
13. rotational velocity 
15. lift curve slope 
17. collective pitch 
19. vert projected area 
21. wing span 
23. wing CDo 
25. horizontal tail area 
27. horizontal tail CL 
29. vertical tail area 
31. vertical tail CL 
33. auxiliary thrust’) 
') 
0. NO CHANGES’ ) 


choice=input (’ 

if isempty (choice), 
choice=0; 

end 

if choices=1, 
cle 


tmp=PA; 
PAzinput (’ Pressure al 
if isempty (PA), 
PA=tmp; 
end 
elseif choice==2, 
clc 
aisp(’ ’) 
temp 
tmp=temp; 
temp=input (’Temperatu 
if isempty (temp), 
temp=tmp; 
end 
elseif choice==3, 
cle 
disp(’ ‘) 
Vinf=Vinf/1.69 
tmp=Vinf ; 


Vinf=input (’Airspeed (knots): ’)*1.69; 


if isempty (Vinf) , 
Vinf=tmp*1.69; 
end 
elseif choice==4, 
ele 


wee’) 


2. 
4. 
6. 
8. 
10. 
12. 
14. 
16. 
18. 
20. 
22. 
24. 
26. 
28. 
30. 
32. 


temperature’ ) 

gross weight’) 

blade radius’ ) 

hinge offset’ ) 

blade twist ’) 

# blade elements’ ) 

# azimuth sectors’ ) 
airfoil’) 

flatplate area’) 

wing area’) 

wing CL’) 

wing efficiency factor’) 
horizontal tail span’) 
horizontal tail CDo’) 
vertical tail span’) 
vertical tail CDo’) 


Input the parameter to change: '); 


titude (ft): 


re (deg F): 
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Gisp(’ ’) 
Gw 


tmp=GW; 
GWsinput (‘Aircraft gross weight (lbs): ’); 
if isempty (GW), 


GWstmp; 
end 
elseif choicess5, 
ele 
disp(’ ‘) 
b 
tmp=b; 


bzinput (‘Number of blades: ‘); 
if isempty (b), 
b=tmp; 
end 
elseif choice=#=:6, 
ele 
Gisp(’ ‘) 
R 


tmp=R; 
R=input (’Blade radius: center of hub to blade tip (ft): ’); 
if isempty (R), 
R=tmp; 
end 
elseif choice==7, 
cle 
disp(’ ‘) 
cblade 
tmp=cblade; 
cblade=input (‘Blade chord (ft): '); 
if isempty (cblade) , 
cblade=tmp; 
end 
elseif choice==8, 
clc 
disp(’ ") 


e 
tmp=e; 
e=input (‘Hinge offset (ft): ‘); 
if isempty(e), 
e=tmp; 
end 
elseif choice==-9, 
cle 
disp(’ ') 
grip 
tmpzgrip; 
gripzinput (’Non-aerodyn inboard portion of blade (ft): ‘); 
if isempty (grip), 
grip=tmp; 
end 
if grip < 1le-10, 
grip-le-10; 
end 
elseif choicez<10, 
cle 
disp(’ ’) 
twistx-twist*S7 .3 
tmp=twist; 
twistsinput (’Blade twist (deg): ‘); 
if isempty (twist), 











twistsabs (tmp) /57.3; 
else 
twisteabs (twist) /57.3; 
end 
elseif choice=<«11, 


ele 
disp (°°) 
wbhlade 
tmp=wblade ; 


wbhladexinput (’Weight of aero portion of blade (lbs): 


if isempty (PA), 
wbhlade=tmp ; 
end 
elseif choicez=#12, 


tmp=nbe ; 
nbesinput (‘Number of blade elements: 
if isempty(nbe), 

nbe=tmp; 
end 

elseif choice==13, 

elec 
disp (’ 
omega 
tmpzomega ; 
omegazinput (‘Rotor rotational velocity (rad/sec) : 
if isempty (omega), 

omega=tmp ; 


elseif choice==14, 
cle 
disp(’ ’) 
naz 
tmp=naz; 
nazzinput (‘Number of azimuth sectors: 
if isempty (naz), 

naz=tmp; 

end 

elseif choice=<=15, 
cle 
disp(’ ‘) 
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") 


oe I 


a 
tmp=a; 
azinput (‘Average lift curve slope (CL vs alpha): 
if isempty (a), 
astmp; 
end 
elseif choice==16, 
cle 
disp(’ ‘) 
airfoil 
tmp=airfoil; 
flag=1 7 
while flag > 0 
airfoilzinput (’Airfoil 
if isempty (airfoil), 
airfoilztmp; 
end 
if airfoil«<1, 
flag=0; 


1. HH-02 2. VR-12 
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1); 


‘); 


>> '); 





elseif airfoil=s2, 

£lag=0; 
else 

disp(’ ’) 

Gisp(’ *** Enter al or 2 ***’) 
end 


end 

elseif choicez=17, 
cle 
disp(' ‘) 
thetaoethetao*57 .3 
tmp=sthetao; 
thetaosinput (‘Collective pitch at .7 r/R (deg): ’)/57.3; 
if isempty (thetao), 

thetao=etmp/57 .3; 

end 

elseif choice==18, 


tmp=Afh ; 
Afhsinput (‘Aircraft equivalent flatplate area (ft*2): ‘'); 
if isempty (Afh), 
Afhz=tmp; 
end 
elseif choice=<=19, 


Afveinput (‘Vertical projected area (ft*2): '); 
if isempty (Afv), 
Afvetmp; 
end 
elseif choice==20, 
cle 
disp(’ ‘) 
Swing 
tmp=Swing ; 
Swing=input (‘Wing area (ft*2): '); 
if isempty (Swing), 
Swingstmp; 
end 
if Swing < le-10, 
Swing=zle-10; 
end 
elseif choicez=«21, 
clc 
disp(’ ’) 
bwing 
tmp=bwing; 
bwing=input ('Wing span (ft) : 
if isempty (bwing) , 
bwingstmp; 
end 
if bwing < le-10, 
bwing=le-10; 
end 
elseif choice=<22, 
cle 
disp(’ ’) 
CLwing 





tmp=CLwing; 
CLwingsinput (‘Expected CL for the wing: ‘); 
if isempty (CLwing) , 
CLwingstap; 
end 
elseif choices«23, 


CDowingsinput (‘Wing profile drag coef (CDo): '); 
if isempty (CDowing), 


CDowing=tmp ; 
end 
elseif choices=«24, 
cle 
disp (‘ *) 
ewing 
tmp=zewing; 


ewingzinput (‘Wing efficiency factor (e): ‘); 
if isempty (ewing) , 
ewingzstmp; 
end 
if ewing < le-10, 
ewingsle-10; 


end 

elseif choice=«25, 
cle 
disp(’ ‘) 
Shoriz 
tmp=Shoriz; 


Shoriz=input (‘Horizontal tail area (ft*2): '); 
if isempty (Shoriz), 
Shoriz=tmp; 
end 
if Shoriz < le-10, 
Shoriz=le-10; 
end 


elseif choice==26, 


cle 

disp(’ ’) 

bhoriz 

tmp=bhoriz; 

bhoriz=input (‘Horizontal tail span (ft): ‘); 

if isempty (bhoriz) , 

bhoriz=tmp; 

end 

if bhoriz < le-10, 
bhoriz=le-10; 

end 


elseif choice==27, 


cle 

disp(’ ’) 
CLhoriz 
tmp=CLhoriz; 


CLhoriz=input (‘Expected CL for the horizontal tail: 


if isempty (CLhoriz) , 
CLhoriz=tmp; 
end 


elseif choice=<28, 


ele 
disp(’ ’) 








'); 








CDohoriz 
tmpeCDohoriz; 
CDohorizsinput (’Horizontal tail profile drag coef (CDo): ‘); 
if isempty (CDohoriz), 
CDohoriz=tmp; 
end 
elseif choices=29, 
cle 
disp(’ ‘) 
Svert 
tmp=Svert ; 
Sverteinput (‘Vertical tail area (ft*2): '); 
if isempty (Svert), 
Svertatmp; 
end 
if Svert < le-10, 
Svertzle-10; 
end 
elseif choice==30, 
cle 
disp(’ ‘) 
bvert 
tmp=bvert ; 
bvert=input (‘Vertical tail span (ft): '); 
if isempty (bvert), 
bvert=tmp; 
end 
if bvert < le-10, 
bvertzle-10; 
end 
elseif choice==31, 


tmp=CLvert ; 
CLvert=input (’Expected CL for the vertical tail: '); 
if isempty (CLvert), 
CLvert=tmp; 
end 
elseif choice==32, 
cle 
disp(’ ’) 
CDovert 
tmp=CDovert ; 
CDovert=input (‘Vertical tail profile drag coef (CDo): '); 
if isempty (CDovert) , 


CDovert=tmp ; 
end 
elseif choice=<=33, 
cle 
disp(’ ‘) 
Taux 
tmp=Taux; 


Tauxsinput (‘Auxiliary thrust (lbs): ’); 
if isempty (Taux), 


Taux=tmp ; 
end 
elseif choicesz=z0, 
elec 
@Gisp(’ ‘) 
disp(’ ’) 
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disp (‘ *** SAVE INSTRUCTIONS ***’) 


disp ( ‘ ’ ) 
disp (‘ A. Save the new data to a specified file name.’) 
Gisp (’ B. Do not use an extension or quotations.‘ ) 
disp (’ C. Use letter/number combinations of 6 characters or less.') 
disp (’ D. The file will be saved with a ".mat" extension.’‘) 
@isp(’ ‘) 
disp (’ ex: dsgn_2’) 
disp(’ ‘) 
disp (’ E. If you made no changes, press < Enter >, the file will’) 
disp (‘ be saved with the original name.’) 
flag=1; 
while flag > 0 
filename0=filenamel ; 
filenamelzinput (’ save file as: ‘','s’); 
if isempty (filenamel), 
filenamel=filenameod; 
end 
clear filenameOd 
if length(filenamel) > 6, 
disp(’ ‘) 
disp (’ use 6 characters or less‘) 
flagz=1; 
else 
flag=0; 
end 
end 
eval([’save ‘',filename1] ) 
check=0; 
else 
disp(’ ') 
disp(’ Enter a displayed number ...press any key to continue’) 
pause 
end 
end 
% ***x If creating a new file: get input for required variables 
% and save under desired file name. Saves to current 
% directory as a .mat file. *** 


elseif answer==2, 
cle 
PA=input (’Pressure altitude (ft): ’); 
while isempty (PA), 
disp(’ ‘) 
disp(’You must enter a numerical value’ ) 
PAzinput (’Pressure altitude (ft): ‘); 
end 
temp=sinput (‘Temperature (deg F): ’); 
while isempty (temp) , 
disp(’ ’) 
disp(’You must enter a numerical value’) 
tempxinput (’Temperature (deg F): ‘); 
end 
Vinf=input (‘Airspeed (knots): ’)¥*1.69; 
while isempty (Vinf), 
disp(’ ’) 
disp(’You must enter a numerical value’) 
vinf=input (‘Airspeed (knots): ’)*1.69; 
end 
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GWainput (‘Aircraft gross weight (lbs): ‘); 
while isempty (GW) , 
disp(’ °) 
disp (’You must enter a numerical value’ ) 
GWeinput (‘Aircraft gross weight (lbs): ‘); 
end 
beinput (‘Mawber of blades: ’); 
while isempty (b), 
disp(’ ‘) 
disp(’You must enter a numerical value’ ) 
bzinput (‘Number of blades: ‘); 
end 
Reinput (’Blade radius: center of hub to blade tip (ft): ‘); 
while isempty (R), 
disp(’ ’) 
disp(’You must enter a numerical value’) 
Rainput (‘Blade radius: center of hub to blade tip (ft): ‘); 
end 
cbhladezinput ('Blade chord (ft): ‘); 
while isempty (cblade) , 
disp(’ ‘) 
disp(’You must enter a numerical value’) 
chladezinput (‘Blade chord (ft): ’); 
end 
exinput (‘Hinge offset (ft): '); 
while isempty(e), 
Gisp(' ‘) 
disp (‘You must enter a numerical value’) 
ezsinput (’Hinge offset (ft): ’); 
end 
grip=input (’Non-aerodynamic inboard portion of blade (ft): '); 
while isempty (grip), 
disp(’ ‘) 
disp (’You must enter a numerical value’) 
grip=input (’Non-aerodynamic inboard portion of blade (ft): '); 
end 
while grip < le-10, 
grip=le-10; 
end 
twist=input (‘Blade twist (deg): ‘)/57.3;,twist=abs (twist) ; 
while isempty (twist) , 
disp(’ ’) 
disp (’You must enter a numerical value’) 
twist=input (’Blade twist (deg): ‘')/57.3;,twist=abs (twist) ; 
end 
wblade=input (‘Weight of aero portion of one blade (lbs): ’); 
while isempty (wblade) , 
disp(’ ’) 
Gisp(’You must enter a numerical value’ ) 
ae de it a of aero portion of one blade (lbs): '); 
en 
nbe=input (‘Number of blade elements: ’); 
while isempty (nbe) , 
disp(’ ’) 
disp(’You must enter a numerical value’) 
nbexinput (‘Number of blade elements: ') ; 
end 
omegazinput (‘Rotor rotational velocity (rad/sec): ‘); 
while isempty (omega) , 
disp(’ ‘) 
disp(’You must enter a numerical value’) 
omega=input (‘Rotor rotational velocity (rad/sec): '); 
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end 
nazzinput (‘Number of azimuth sectors: '); 
while isempty (naz), 
disp (’ *) 
disp(’You must enter a numerical value’) 
nazesinput (’Number of azimuth sectors: '); 
end 
azinput (‘Lift curve slope of rotor airfoil (CL vs alpha): ‘); 
while isempty (a), 
disp(’ ’) 
disp (’You must enter a numerical value’) 
azinput (‘Lift curve slope of rotor airfoil (CL vs alpha): ’); 
end 
flag=1; 
while flag > 0 
airfoilzinput (‘Airfoil 1. HH-02 2. VR-12 >> '); 
while isempty (airfoil), 
disp(’ ’) 
disp (‘You must enter a numerical value’) 
airfoilzinput (’Airfoil 1. HH-02 2. VR-12 >> '); 
end 
if airfoil~=1 & airfoil~=2, 
disp(’ ‘) 
Gisp(’ *** Enter al or 2 ***’) 
else 
flag=0; 
end 
end 
thetaozinput (‘Collective pitch at .7 r/R (deg): ')/57.3; 
while isempty (thetao), 
disp(’ ’) 
disp(’You must enter a numerical value’) 
thetao=input (‘Collective pitch at .7 r/R (deg): ’)/57.3; 
end 
Afhzinput (/Aircraft equivalent flatplate area (ft*2): ’); 
while isempty (Afh), 
disp(’ ’) 
disp(’You must enter a numerical value’) 
Afhsinput (’Aircraft equivalent flatplate area (ft*2): ’); 
end 
Afvzinput (’Vertical projected area (ft*2): '); 
while isempty (Afv) , 
disp(’ ’) 
disp(’You must enter a numerical value’) 
Afvzinput (’Vertical projected area (ft*2): '); 
end 
Swing=input (’Wing area, 0 if no wing (ft*2): ’); 
while isempty (Swing) , 
@isp(’ ’) 
disp(’You must enter a numerical value’) 
Swing=zinput (’Wing area, 0 if no wing (ft*2): ’); 
end 
if Swing~=0, 
bwing=input (‘Wing span (ft): ‘); 
while isempty (bwing), 
disp(’ ’) 
disp (’You must enter a numerical value’) 
bwing=input (’Wing span (ft): ’); 
end 
if bwing < le-10, 
bwing=1le-10; 
end 





CLwingsinput (‘Expected CL for the wing: ‘); 
while isempty (CLwing) , 
disp(’ *) 
Gisp (‘You must enter a numerical value’) 
CLwingswinput (‘Expected CL for the wing: ‘); 
end 
CDowingzinput (‘Wing profile drag coef (CDo): ’); 
while isempty (CDowing) , 
disp(’ ‘) 
disp(’You must enter a numerical value’) 
CDowingszinput (‘Wing profile drag coef (CDo): ‘); 
end 
ewingzinput (‘Wing efficiency factor (e): ’); 
while isempty (ewing) , 
disp(’ ‘) 
disp(’You must enter a numerical value’) 
ewingzinput (‘Wing efficiency factor (e): ‘'); 
end 
if ewing < le-10, 
ewing=1le-10; 
end 
else 
Swing=le-10; 
bwing=le-10; 
CLwinge0; 
CDowing=0; 
ewing=1le-10; 
end 
Shoriz=input (’Horizontal tail area, 0 if none (ft*2): ‘); 
while isempty (Shoriz), 
disp(’ ’) 
disp(’You must enter a numerical value’) 
Shorizzinput (‘Horizontal tail area, 0 if none (ft*2): '); 
end 
if Shoriz~=0, 
bhoriz=input (‘Horizontal tail span (ft): ’): 
while isempty (bhoriz), 
disp(’ ') 
disp(’You must enter a numerical value’) 
bhoriz=sinput (‘Horizontal tail span (ft): ’); 
end 
if bhoriz «< 1le-10, 
bhoriz=le-10; 
end 
CLhoriz=input (‘Expected CL for the horizontal tail: '); 
while isempty (CLhoriz) , 
disp(’ ‘) 
disp(’You must enter a numerical value’) 
CLhoriz=input (’Expected CL for the horizontal tail: ‘); 
end 
CDohorizzinput (‘Horizontal tail profile drag coef (CDo): ‘); 
while isempty (CDohoriz), 
Gisp(’ ’) 
disp(’You must enter a numerical value’ ) 
CDohoriz=input (‘Horizontal tail profile drag coef (CDo): '); 
end 
else 
Shoriz=le-10; 
bhoriz=#le-10; 
CLhoriz=0; 
CDohoriz=0; 
end 
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Svertzinput (‘Vertical tail area, 0 if none (ft*2): ‘); 
while isempty (Svert) , 
diap(’ ') 
disp(’You must enter a numerical value’) 
Svertszinput (’Vertical tail area, 0 if none (ft*2): ‘); 
end 
if Svert~=0, 
bvertzinput (‘Vertical tail span (ft): '); 
while isempty (bvert), 
disp(’ ’) 
disp(’You must enter a numerical value’) 
bvert=input ('Vertical tail span (ft): ’); 
end 
if bvert < le-10, 
bvertsle-10; 
end 
CLvertzinput (‘Expected CL for the vertical tail: ’); 
while isempty (CLvert) , 
disp(’ ‘) 
disp(’You must enter a numerical value’) 
CLvert=input (‘Expected CL for the vertical tail: '); 
end 
CDovert=input (’Vertical tail profile drag coef (CDo): ‘); 
while isempty (CDovert) , 
disp(’ ‘) 
Gisp(’You must enter a numerical value’) 
CDovert=input (’Vertical tail profile drag coef (CDo): ‘); 
end 
else 
Svert=le-10; 
bvert=le-10; 
CLvert=0; 
CDovert=0; 
end 
Tauxsinput (‘Auxiliary thrust (lbs): ’); 
while isempty (Taux) , 
disp(’ ’) 
disp(‘’You must enter a numerical value’ ) 
Taux=input (‘Auxiliary thrust (lbs): ’); 


end 
elc 
disp (' ") 
disp(’ ’) 
disp (’ *#** SAVE INSTRUCTIONS ***’) 
disp(’ ‘) 
disp (’ A. Save the data to a specified file name.’) 
disp (’ B. Do not use an extension or quotations.’ ) 
disp (’ C. Use letter/number combinations of 6 characters or less.’) 
disp (’ D. The file will be saved with a ".mat" extension.’ ) 
@Gisp(’ ') 
disp (’ ex: dsgn_1’) 
Gisp(’ ’) 
disp (’ E. If you do not enter a name, the default is "dsgn_1" ’) 
flag=1; 
while flag > 0 
filenamel=input (’ save file as: ','8s'); 


if isempty (filenamel), 
filenamel='dsgn_1’; 

end 

if length(filenamel) > 6, 
disp(’ °) 
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disp (’ use 6 characters or less’) 
flags1; 
else 
flag=0; 
end 
end 
eval((’save ‘,filename1) ) 
check=0; 
else 
disp(’ ’) 
Gisp (’ Enter a 1 or 2’) 


*** DATA ENTRY COMPLETE ***’) 
pause (3) 


check2=1; 
while check2 > 0 


cle 
disp(’ ‘) 
disp (’ *** EXECUTION MENU ***’) 
disp(’ ') 
disp (’ . Rotor Performance Analysis’) 
disp (' . Stability and Control Analysis’) 
disp (’ . Rotor Dynamics Analysis’ ) 
disp (’ . Change data’) 
disp (’ - Quit’) 
check3=1; 
while check3 > 0 
answer=input (’ Enter al, 2, 3, 4, or 5 
save temp checkl check3 filenamel 


elseif answer 
stab 
load temp 
check3=0; 


elseif answer 3, 

insert name blade dynamics M-file here 

load temp 

check3=0; 
disp(’ ‘) 
disp (’ Blade dynamics module is currently under development’ ) 
disp(’ ') 
disp (’ press any key to continue...’) 
pause 


elseif answer 4, 
cle 
check2=0; 
check3=0; 








elseif answer == 5, 
disp(’ ‘) 
check1s0; 
check2=0; 
check3=0; 

else 
disp(’ ‘) 

end 

end 
end 
end 





APPENDIX B: TRIM.M 


*% Rotor trim subroutine 


disp(’ ’) 
disp (' *** EXECUTING ROTOR TRIM ROUTINE ***’) 


% *** calculation of required parameters *** 
rhoz .002377* (- .000031*PA+(-.002*temp+1.118) ) ; 
% **« first guess at rotor profile drag ( H force) *** 


if Vint < 16.9, 
Drotor=0; 
else 
Drotor=Vinf* (rho/ .002377) ; 
end 


q=0.5*rho*vinf*2; 

Adisk=pi*R‘*2; 

Vtip=omega*R; 

Dfuse=q*Afh; 
CDwing=CDowing+ (CLwing*2/ (ewing*pi* (bwing*2/Swing) )) ; 
CDhori.z=CDohoriz+ (CLhoriz*2/(.8*pi* (bhoriz*2/Shoriz) )); 
CDvert=CDovert+ (CLvert*2/(.8*pi* (bvert*2/Svert) )); 
Dwing=q*CDwing*Swing; 

Dhoriz=q*CDhoriz*Shoriz; 

Dvert=q*CDvert*Svert ; 

Dftotal= (Dfuse+Dwing+Dhoriz+Dvert) -Taux; 
Lwing=q*CLwing*Swing; 

Lhoriz=q*CLhoriz*Shoriz; 

Lvert=q*CLvert *Svert ; 

Lftotal=Lwing+Lhoriz; 
alphaT=atan2 ( (Dftotal+Drotor) , (GW-Lftotal) ) ; 
musVinf*cos (alphaT) /Vtip; 


% *** thrust calculation *** 


if Vinf < 16.9, 
T= (1+ (0.3*A£v/Adisk) ) *GW; 
CT=T/ (Adisk*rho*Vtip*2) ; 
else 
T= (GW-Lftotal) /cos (alphatT) ; 
CT=T/ (Adisk*rho*Vtip“*2) ; 
end 


* *** setup blade radius elements, azimuth elements, 


% induced velocity distributions, and determination 
% ef coning angle and tip loss parameter *** 

Bel- (sqrt (2*CT) /b) ; 

Reff=<B*R; 

Rbar=Reff-e; 


dr= (Reff-grip) /nbe; 
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regrip:dr:Reff-dr;,r=sr+dr/2; 

xT1<0.7;,% *** first guess at rT *** 

RbarT=rT1*Rbar; 

mblade=wblade/32.17; 

betaozasin ( (T/b*RbarT- (.5* (R-e) +e) *wblade) /((.5* (R-e) +e) *2*omega*2*mblad 
e)); 

betatetwist* (0.7-(r/R)); 

psi=0:360/naz:360-360/naz; ,psixpsi’ /57.3; 


% «** induced velocity determination *** 


if Vinf < 16.9, 

A=4*pi; 

Bvz (b/2) *omega*a*cblade; 

Tv=0; 

delT=T-Tv; 

while abs(delT) > .01*T 
thetavstwist* (0.7- (r/R) )+thetao; 
C= (-b/2) *cblade*omega*2*r*a.*thetav; 
vis (-Bv+sqrt (Bv*2- (4*A*C) )) /(2*A) ; 


€aTv= (b/2) *rho* ( (omega*r) .*2) *a.* (thetav- (vi. / (omega*r) )) *cblade*dr; 
Tv=sum (daTv) ; 
GelT=T-Tv; 
if delT < 0, 
thetao=thetao-0.5*thetao*abs (delT/T) ; 
else 
thetao=thetao+0.5*thetao*abs (delT/T) ; 
end 
end 
else 
lamdaT=[1 -2*mu*sin (alphaT) 
mu*2* (sin (alphaT) *2+1) -2*mu*3*sin(alphaT) mu*4*sin(alphaT) *2-CT*2/4] ; 
lamdaT=max (real (roots (lamdaT) ) ) ; 
vislamdaT*Vtip-Vinf*sin(alphaT) ; 
vizvi*ones (r) ; 
end 


% *** first guess at theta *** 


thetalc=0 .035* ((0.0006e-3*Vinf*2+0 .244e-3*Vinf) /0.105) ; 
thetals=-0.087* ((0.0006e-3*Vinf*2+0.244e-3*Vinf) /0.105) ; 
theta=thetao+thetalic.*cos (psi) +thetals.*sin (psi) ; 


% *** rotor trimming routine *** 


disp (’ *** TRIMMING COLLECTIVE ***’) 


k=1; 
error0=(T*.02) +41; 


while abs(error0) > T*.02 
Tpsi=zeros (psi) ; 
threalc 
error0=T- (mean (Tpsi) *b) ; 
if errord «< -T*.02, 
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thetaosthetao-0.35*thetao*abs (1.5*error0/T) * (1-mu) ; 


elseif error0 > T*.02, 


thetaozethetao+0.35*thetao*abs (1 .5*error0/T) * (1-mu) ; 


end 
thetaz=thetao+thetalc. *cos (psi) +thetals.*sin (psi) ; 
if k > 1, 


if abs(error0) > abs(errorl), 

cle 

aisp(’ ‘) 

disp(’ ‘) 

disp(’This configuration will not trim’) 
disp(’Try a lower airspeed or a new design’) 


disp(’ ‘) 
disp(’ Program execution terminated’ ) 
disp(’ ‘) 
error ('*** END OF PROGRAM ***’) 
end 
errorlzerror0; 
k=xk+1 3 


") 
wwe TRIMMING CYCLIC ***’) 


t0=clock; 


error0=( ((T/b) *rT1* (R-grip) )*.04) 41; 


while errorO > ((T/b) *rT1*(R-grip))*.04 


time=etime (clock, t0) ; 
if time > 15, 


disp(’ ‘) 
disp(’still trimming ...’) 
t0=clock; 


Mpsi(:,k) =zeros (psi) ; 
tmcalc 

theta= (theta theta(:,k)]; 
Mpsi=(Mpsi Mpsi(:,k)]; 


*** calculation of initial dthetadM *** 


if k < 2, 

theta(:,k+1) =theta(:,k)+0.25/57.3; 

Mpsi(:,k+1) =zeros (psi) ; 

k=k+1; 

tmealc 

k=k-1; 

athetadM= (theta(:,k+1) -theta(:,k)) ./(Mpsi(:,k+1) -Mpsi(:,k)); 
end 


*** calculation of M first harmonic parameters *** 


Mic=2*sum(Mpsi (:,k) .*cos (psi) ) /naz; 
Mls=2*sum(Mpsi(:,k) .*sin(psi)) /naz; 


*“** yemoval of first harmonic terms from Mpsi *** 
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BO 


Mpsi(:,k+1) sMpsi (:,k) -Mic.*cos (psi) -Mls.*sin (psi) ; 
delM=Mpsi (:,k+1)-Mpsi(:,k); 
error0smax (de1M) -min (delM) ; 


if k > 1, 
if error0O > errorl, 
cle 
disp(’ ’) 
disp(‘ ’) 


disp(’This configuration will not trim’) 
disp(’Try a lower airspeed or a new design’) 
Gisp(’ ’) 
G@isp(’ Program execution terminated’ ) 
disp(’ ’) 
error (‘*** END OF PROGRAM ***’) 
end 
end 
errorl=serroro; 


* *** calculation of new theta *** 


GelM=0 .5* (1-mu) *delM; 
theta(:,k+1) =theta(:,k)+(dthetadM.*delM) ; 


if errorO <= ((T/b) *rT1* (R-grip))*.04, 


theta.. ‘sum(theta(:,k) .*cos (psi) ) /naz; 
thetals=<*sum(theta(:,k) .*sin(psi) ) /naz; 
else 


thetalc=2*sum(theta(:,k+1) .*cos (psi) ) /naz; 
thetals=2*sum(theta(:,k+1) .*sin(psi) ) /naz; 
end 


theta(:,k+1) =thetao+thetalc.*cos (psi) +thetals.*sin (psi) ; 
% *** calculation of new dthetadM *** 


theta= [theta theta(:,k+1)]; 
Mpsi=([Mpsi Mpsi(:,k+1)); 
theta(:,k+2) =theta(:,k)+0.25/57.3; 
Mpsi (:,k+2) =zeros (Mpsi(:,k+1)); 
k=k+2 ? 
tmcale 
k=k-2; 
dthetadM= (theta(:,k+2)-theta(:,k)) ./(Mpsi(:,k+2) -Mpsi(:,k)); 
k=k+1 , 
end 


disp(’ ‘) 
disp (’ *** ADJUSTING COLLECTIVE ***’) 
disp(’ ‘) 


theta=theta(:,k); 
kel; 
error0=(T*.01) +1; 


while abs(error0) > T*.01 
Tpsi=zeros (psi) ; 
threalc 
error0=T- (mean (Tpsi) *b) ; 
if error0 < -T*.01, 
thetao=thetao-0.25*thetao*abs (1.25*error0/T) * (1-mvu) ; 
elseif error0O > T*.01, 
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thetaosthetao+0 .25*thetao*abs (1.25*error0/T) * (1-mu) ; 
end 
thetaethetao+thetalc.*cos (psi) +thetals.*sin (psi) ; 
if k > 1, 
if abs(error0) > abs(errorl), 


disp(’ ‘) 
disp (’This configuration will not trim’) 
disp(‘Try a lower airspeed or a new design’) 
Gisp(’ ‘) 
disp (’ Program execution terminated’) 
disp(’ ‘) 
error (‘*** END OF PROGRAM ***’) 
end 

end 

errorlzerroro0; 

kzk+1 3 

end 


% *** calculating drag moments *** 


DMpsi=zeros (psi) ; 
amcalc 


% *** calculating rotor H force *** 


if Vinf < 16.9, 
Hrotors0; 
aT= [dT ddT] ; 
aD=[dD ddD] ; 
else 
aT=[AT ddT] ; 
aD=(dD ddbD]; 
for isl:length (r) +1, 
Hlc(i) =2*sum (dT(:,1i) .*cos (psi) ) /naz; 
H1s (i) =2*sum(dD(:,i) .*sin(psi)) /naz; 
end 


Hrotors (( (b*cos (alphaT) /2) * (sum(H1s) -sin (betao) *sum (Hic) ))+Drotor) /2; 
end 
% *** calculating new rT *** 


rT2= (( (mean (Mpsi (:, length (Mpsi (1, :))-1)) /mean(Tpsi) ) /R) +rT1) /2; 


% *** check rotor drag and rT, retrim rotor if required *** 
while abs(Drotor-Hrotor) > 0.2*Hrotor | abs(rT1l-rT2) > 0.015*rT1l 


if abs (Drotor-Hrotor) > 0.2*Hrotor, 
Gisp(’ ’) 
disp (' *** ADJUSTING ROTOR DRAG ***’) 
end 
Drotor=Hrotor; 
if abs(rT1l-rT2) > 0.015*rT1, 
disp(’ ’) 
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disp (’ *** ADJUSTING MEAN THRUST LOCATION ***’) 
end 


disp(’ ‘) 
disp (’ *** RETRIMMING ROTOR ***’) 
disp(’ ‘) 
OT=#dT(:,1: 
@DadD(:,1 


% *** recalculating parameters *** 


alphaTzatan ( (Dftotal+Drotor) / (GW-Lftotal) ) ; 
musVinf*cos (alphaT) /Vtip; 


if Vinf >= 16.9, 
T= (GW-Lftotal) /cos (alpha?) ; 
CT=T/ (Adisk*rho*Vtip~2) ; 
lamdaTz[1 -2*mu*sin (alphaT) 

mu*2* (gin (alphaT) *2+1) -2*mu*3*sin(alphaT) mu*4*sin(alphaT) *2-CT*2/4) ; 

lamdaT=max (real (roots (lamdaT) ) ) ; 
vizlamdaT*Vtip-Vinf*sin (alphaT) ; 
vizvi*ones (r) ; 

end 


B=1- (sqrt (2*CT) /b) ; 
Reff=B*R; 

Rbar=Reff-e; 

dr= (Reff-grip) /nbe; 
r=grip:dr:Reff-dr;,r=r+dr/2; 
RbarT=rT2*Rbar ; 


betaozasin ( (T/b*RbarT- (.5* (R-e) +e) *wblade) /((.5* (R-e) +e) “2*omega*2*mblad 
e)); 


% *** trimming collective *** 


tO=clock; 
kel; 
error0=(T*.02) +41; 


while abs(error0) > T*.02 
Tpsi=zeros (psi} ; 
threalec 
error0=T- (mean (Tpsi) *b) ; 
if error0O < -T*.02, 
thetao=thetao-0.35*thetao*abs (1.5*error0/T) * (1-mvu) ; 
elseif errorO > T*.02, 
thetao=thetao+0.35*thetao*abs (1.5*error0/T) * (1-mu) ; 
end 
theta=thetao+thetalc.*cos (psi) +thetals.*sin (psi) ; 
if k > 1, 
if abs(error0) > abs(errorl), 
eclc 
disp(’ ‘) 
Gisp(’ ’) 
disp (‘This configuration will not trim’) 
disp(’/Try a lower airspeed or a new design’) 
disp(’ ’) 
disp (‘Program execution terminated’ ) 
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disp(‘ ‘) 
error (’*** END OF PROGRAM ***’) 
end 
end 
erroriszserroro; 
kwxk+1; 
end 


% *** trimming cyclic *** 


k=1; 
errors (((T/b) *rT2* (R-grip) ) *.04) +1; 


while errorO > ((T/b) *rT2* (R-grip) )*.04 


timezetime (clock, t0) ; 

if time > 15, 
disp(’still trimming ...’) 
disp(’ ‘) 
tO#clock; 

end 


Mpsi(:,k) =zeros (psi) ; 
tmcalc 

theta=(theta theta(:,k)]; 
Mpsi=([Mpsi Mpsi(:,k)]; 


*** calculation of initial dthetadM *** 


if k < 2, 

theta (:,k+1) =theta(:,k)+0.25/57.3; 

Mpsi (:,k+1) =zeros (psi) ; 

kz=k+1; 

tmcalc 

k=k-1; 

athetadM= (theta(:,k+1)-theta(:,k)) ./(Mpsi(:,k+1) -Mpsi(:,k)); 
end 


*** calculation of M first harmonic parameters *** 


Mlc=2*sum (Mpsi(:,k) .*cos (psi) ) /naz; 
Mls=2*sum (Mpsi (:,k) .*sin (psi) ) /naz; 


*** removal of first harmonic terms from Mpsi *** 


Mpsi (:,k+1) =Mpsi (:,k) -Mic.*cos (psi) -Mis.*sin (psi) ; 
delM=Mpsi (:,k+1) -Mpsi(:,k); 
error0smax (de1M) -min (de1M) ; 
if k > 1, 
if errorO > errorl, 
elec 
disp (’ ") 
disp(’ ‘) 
disp(’This configuration will not trim’) 
disp(’Try a lower airspeed or a new design’) 
Gisp(’ ‘) 
disp(’ Program execution terminated’ ) 
disp(’ ‘) 
error (‘'*** END OF PROGRAM ***') 
end 





errorliserror0o; 
*** calculation of new theta *** 


delM=0 .5* (1-mu) *delM; 

theta(:,k+1) «theta(:,k)+(dthetadM.*delM) ; 

if error0O <= ((T/b) *rT2* (R-grip))*.04, 
thetalc=2*sum(theta(:,k) .*cos (psi) ) /naz; 
thetals=2*sum(theta(:,k) .*sin(psi) ) /naz; 

else 
thetalc=2*sum(theta(:,k+1) .*cos (psi) ) /naz; 
thetals«2* sum (theta(:,k+1) .*sin(psi) ) /naz; 

end 

theta(:,k+1) sthetao+thetalc.*cos (psi) +thetals.*sin (psi) ; 


*** calculation of new dthetadM *** 


theta=[theta theta(:,k+1)]; 
Mpsi=[Mpsi Mpsi(:,k+1))]; 
theta(:,k+2) =theta(:,k)+0.25/57.3; 
Mpsi (:,k+2) =zeros (Mpsi(:,k+1)); 
k=k+2 ; 
tmcalc 
kzk-2 : 
dathetadM= (theta(:,k+2) -theta(:,k)) ./(Mpsi(:,k+2) -Mpsi(:,k)); 
k=k+1 , 
end 


*** retrimming collective *** 


thetastheta(:,k) ; 
k=1; 
error0=(T*.01) 41; 


while abs(error0) > T*.01 
Tpsi=zeros (psi) ; 
thrcealc 
error0=T- (mean (Tpsi) *b) ; 
if errorO < -T*.01, 
thetaosthetao-0.25*thetao*abs (1.25*error0/T) * (1-mu) ; 
elseif error0O > T*.01, 
thetaosthetao+0.25*thetao*abs (1.25*error0/T) * (1-mu) ; 
end 
theta=thetao+thetalc.*cos (psi) +thetals.*sin (psi) ; 
if k > 1, 
if abs(error0) > abs(errorl), 


cle 
disp(’ ‘) 
disp(’ ‘) 


disp (‘This configuration will not trim’) 
disp(’Try a lower airspeed or a new design’) 


disp(’ ’) 
disp(’ Program execution terminated’ ) 
disp(’ ’) 
error (’*** END OF PROGRAM ***’) 
end 

end 

errorizerror0o; 

kzk+1 3 

end 
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% *** recalculating rotor H force *** 


if Vinf < 16.9, 
Hrotors0; 
@T=(@T daT) ; 
adD=(dD ddDj; 
else 
aT= (AT adaT) : 
dD=(aD dadD] ; 
for izl:length(r)+1, 
Hlc (i) =2*sum(daT(:,1) .*cos (psi) ) /naz; 
Fe ae aan Cea) oem EOE t)) (BRET 
en 


Hrotors ( ( (b*cos (alphaT) /2) * (sum (H1s) - sin (betao) *sum(H1ic) ))+Drotor) /2; 
end 


% *** recalculating rT *** 


rTl=rT2; 
miaek pean pet te; tengee (pet py ket idmeen tell yin) xe et he: 
en 


* *** recalculating drag moments *** 


dT=dT(:,1:nbe) ; 
@D=dD(:,1:nbe) ; 
DMpsi=zeros (psi) ; 
dmcalc 

aT=[(daT ddT] ; 
@D=[daD ddbD] ; 

cle 


disp (’ *** ROTOR TRIMMED ***’) 
pause (3) 
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APPENDIX C: PERF.M 


% *** Performance output subroutine *** 
% *** Clearing all variables in the Matlab environment *** 
clear 


% *** Loading input and check parameters *** 


load temp 

eval((’load ',filename1) ) 

cle 

disp(' ‘) 

disp (’ **w ROTOR PERFORMANCE ROUTINE ***’) 
disp (’ waeeeknn!) 

disp(’ ‘) 

trim 


% *** Calculation of output parameters *** 


Protor=mean (DMpsi) *b*omega/550; 
Qrotor=mean (DMpsi) *b; 
solidity=b*cblade/ (pi*R) ; 
CQ=Qrotor/ (Adisk*rho*Vtip*2*R) ; 
CH=Hrotor/ (Adisk*rho*Vtip“2) ; 
CT_sig=CT/solidity; 
CQ_sig=CQ/solidity; 
CH_sig=CH/solidity; 
Machtips= (Vtip*cos (alphaT) +Vinf) / (49 .05*sqrt (temp+460) ) ; 
if Vint < 16.9, 
DL=T/ (pi*R*2) ; 
FM= (T* sqrt (DL/ (2*xrho) )) /(550*Protor) ; 
else 
DL=0; 
FM=0; 
end 


ele 
disp(’ ') 
disp (’ *** PERFORMANCE CALCULATIONS COMPLETE ***’) 
Gisp(’ ’) 
disp (’ Do you want the performance results displayed on 
screen?’ ) 
flag=1; 
while flag > 0 

answerzinput (’ 1. yes 2. no >> '); 

if answer == 2, 

flag=0; 
elseif answer == 1, 


% *** output to screen *** 
elc 


disp(’ ’) 
disp (’ *** INPUT DATA ***’) 
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eval ([‘disp(’’ ',filenamel,‘'’)')) 
disp(' ‘) 
fprint£ (’ Forward velocity = %6.0f kts\n’,Vinf/1.69) 
fprintf (’ Temperature = %6.0f degs F\n’,temp) 
fprint£ (’ Pressure altitude = %6.0f ft\n’, PA) 
fprint®£ (’ Gross weight = %6.0£ lbs\n’ ,GwW) 
fprint£ (’ Number of blades = %6.0f \n’,b) 
fprintf (’ Rotor radius = %6.2£ ft\n’,R) 
fprintf (’ Blade chord = %6.2f£ ft\n’,cblade) 
fprint£ (’ Blade twist = %6.2f degs\n’, -1*twist*57.3) 
if airfoil==1, 

disp (’ Blade airfoil = HH-02’) 
else 

disp (’ Blade airfoil = VR-12’) 
end 
fprintf (’ Blade lift curve slope = %6.2f \n’,a) 
fprintf (’ Blade weight = %6.2f lbs\n’,wblade) 
fprint£ (’ Rotational velocity = %6.2f rads/sec\n’ ,omega) 
fprint£ (’ Blade grip length = %6.2f ft\n’',grip) 
fprinté£ (’ Hinge offset = %6.2f ft\n’,e) 
disp(’ ‘) 
disp(’ ‘) 
disp(’press any key to continue...’) 
Pause 
ele 
disp(’ ‘) 
disp (’ *** INPUT DATA CONTINUED ***’) 
eval ([’disp(’’ ‘, filenamel,’’’)’)) 
disp(’ ‘) 
fprintf (’ Equivalent flat plate area = %6.2f ft*2\n’ ,Afh) 
fprintf (' Vertical projected area = %6.2£ f£t*2\n’ ,Afv) 
fprintf (’ Wing area = %6.2f ft*2\n’,Swing) 
fprintf (’ Wing span = %6.2f ft\n’,bwing) 
fprintf (‘ Wing CL = %6.2£ \n’,CLwing) 
fprintf (’ Wing CDo = %6.4£ \n’,CDowing) 
fprintf (’ Wing efficiency factor = %6.2f \n’,ewing) 
fprintf (’ Horizontal tail area = %6.2f ft*2\n’,Shoriz) 
fprinté£ (' Horizontal tail span = %6.2f ft\n’,bhoriz) 
fprintf (’ Horizontal tail CL = %6.2f£ \n’,CLhoriz) 
fprintf (’ Horizontal tail CDo = %6.4f \n’,CDohoriz) 
fprintfé (’ Vertical tail area = %6.2f ft*2\n’,Svert) 
fprintf (’ Vertical tail span = %6.2f ft\n’,bvert) 
fprint£ (’ Vertical tail CL = %6.2£ \n’,CLvert) 
fprintf (' Vertical tail CDo = %6.4f \n’,CDovert) 
fprintf (’ Auxiliary thrust = %6.0f lbs\n’ , Taux) 
disp(’ ‘) 
disp(’ ') 
disp (‘press any key to continue...’) 
pause 
cle 
disp(’ ‘) 
disp (’ *** CALCULATED DATA ***’) 
eval ([’disp(’’ ‘, filenamel,’’’)’)) 
disp(’ ') 
fprintf (’ Fuselage drag = *6.0f lbs\n’ ,Dfuse) 
fprintf (’ Rotor drag = %6.0f lbs\n’ ,Hrotor) 
fprintf (’ Wing lift = %6.0£ lbs\n’,Lwing) 
fprintf (’ Wing drag = %6.0£ 1lbs\n’ ,Dwing) 
fprinté£ (' Horizontal tail lift = %6.0f 1lbs\n’,Lhoriz) 
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%6.0£ lbs\n’ , Dhoriz) 

%6.0£ lbs\n’,Lvert) 

%6.0f lbs\n’ ,Dvert) 

%6.2£ degs\n’ ,alphaT*57.3) 
%6.2£ degs\n’' ,betao*57.3) 
%6.2f \n’,rT2) 

%6.2£ degs\n’ , thetao*57. 3) 
%6.2£ \n’,thetalc*57.3) 
%6.2f \n’,thetals*57.3) 


',filenamei1,’’’)’)) 


%6.3f£ \n’,solidity) 
%6.2£ lbs/ft*2\n’ ,DL) 
%6.2£ \n’,FM) 

%6.3f \n’,CT_sig) 
%6.4£ \n’,CQ sig) 
%6.4f \n’,CH_sig) 
%6.3£ \n’,Machtip) 
%6.3£ \n’,mu) 

%6.0£ lbs\n’,T) 

%6.0f h.p.\n’, Protor) 
%6.0£ ft-1lbs\n’ ,Qrotor) 


fprintf (’ Horizontal tail drag = 
fprintf (’ Vertical tail side force = 
fprintf (' Vertical tail drag = 
fprintf (’ Tip path angle = 
fprinté (‘ Rotor coning angle = 
fprintf (’Location of mean thrust (r/R) = 
fprintt (’ Collective pitch at .7 r/R = 
fprintf(’ 1st lat cyclic term-Al (deg) = 
fprintf(’1lst long cyclic term-Bl (deg) = 
disp(’ ‘) 
disp(’ ’) 
disp(’press any key to continue...’) 
pause 
cle 
dGisp(’ ‘) 
disp (' *** CALCULATED DATA CONTINUED ***’) 
eval ([’disp(‘’ 
disp(‘ ’) 
fprintf (’ solidity (sigma) = 
fprintf (’ Disk loading = 
fprintf (’ Figure of Meri‘ = 
fprintf (’ CT/sigma = 
fprintf (’ CQ/sigma = 
fprinté£ (’ CH/sigma = 
fprintf (’ Tip mach of the adv. blade = 
fprintf£ (’ Advance ratio = 
fprintf(’ Rotor thrust required (TPP) = 
fprintf (’ Rotor power required = 
fprint€ (’ Rotor torque = 
aisp(’ ’) 
disp(’ ’) 
disp(’press any key to continue...’) 
pause 
elec 

flag=0; 

else 
disp(’' ‘’) 
disp (’ Enter a 1 or 2’) 
end 

end 


% *** output to disk (text file) *** 


disp(’ ‘) 
disp (’ 

pause (2) 
diary on 


diary off 


*** Saving data to disk ***’) 


delete diary 


diary on 
disp (’ 


eval ([’disp(’’ 


Gisp(’ ‘) 
fprintf (’ 
fprinté (’ 
fprintf (’ 
fprinté£ (’ 
fprinté£ (’ 
fprintf (’ 


keke RESULTS wen!) 


Forward velocity 
Temperature 
Pressure altitude 
Gross weight 
Number of blades 
Rotor radius 
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‘, filenamel1,’’’)’)) 


%6.0f 
Of 
Of 
Of 
%6.0f 
.2f 


kts\n’, Vinf/1.69) 
degs F\n’,temp) 
ft\n’ , PA) 
lbs\n’ , GW) 

\n’ ,b) 

ft\n’ ,R) 





fprintt (’ Blade chord %6.2f ft\n’,cblade) 


%6.2£ \n’,thetalc*57.3) 
$6.2£ \n’,thetals*57.3) 


fprintf(’ lst lat cyclic term-Al (deg) 
fprintf(’1lst long cyclic term-Bl1 (deg) 


fprintf (’ Blade twist = %6.2f£ degs\n’,-1*twist*57.3) 
fprintf (’ Blade lift curve slope = %6.2f \n’,a) 
fprintf (’ Blade weight = %6.2f lbs\n’ , wolade) 
fprint£ (’ Rotational velocity = %6.2f rads/sec\n’‘ , omega) 
fprintf (‘ Blade grip length = %6.2f ft\n’,grip) 
fprintf (‘ Hinge offset = %6.2f ft\n’,e) 
fprintf (’ Equivalent flat plate area = %6.2f ft*2\n’ ,Afh) 
fprinté£ (’ Vertical projected area = %6.2f£ ft*2\n’ ,Afv) 
fprintf (’ Wing area = %6.2f ft*2\n’,Swing) 
fprintf (’ Wing span = %6.2f ft\n’,bwing) 
fprintf (’ Wing CL = %6.2f£ \n’,CLwing) 
fprintf (' Wing CDo = %6.4f£ \n’,CDowing) 
fprintf (’ Wing efficiency factor = %6.2f \n’,ewing) 
fprintf (’ Horizontal tail area = %6.2f £t*2\n’,Shoriz) 
fprintf (’ Horizontal tail span = %6.2f ft\n’,bhoriz) 
fprintf (’ Horizontal tail CL = %6.2f \n’,CLhoriz) 
fprintf (’ Horizontal tail CDo = %6.4f \n’,CDohoriz) 
fprintf (’ Vertical tail area = %6.2f ft*2\n’,Svert) 
fprintf (’ Vertical tail span = %6.2£ ft\n’,bvert) 
fprintf (’ Vertical taii CL = %6.2f \n’,CLvert) 
fprintf (’ Vertical tail CDo = %6.4£ \n’,CDovert) 
fprintf (' Fuselage drag = %6.0f 1lbs\n’ ,Dfuse) 
fprintf£ (’ Rotor drag = %6.0f lbs\n’ ,Hrotor) 
fprintf (’ Wing lift = %6.0£ lbs\n’, Lwing) 
fprintf (’ Wing drag = %6.0f lbs\n’ ,Dwing) 
fprintf Horizontal tail lift = %6.0£ lbs\n’ , Lhoriz) 
fprintf , Horizontal tail drag = %6.0f lbs\n’ ,Dhoriz) 
fprintf (’ Vertical tail side force = %6.0f£ lbs\n’,Lvert) 
fprintf£ (’ Vertical tail drag = %6.0f lbs\n’ ,Dvert) 
fprintf (’ Auxiliary thrust = %6.0f lbs\n’ ,Taux) 
fprintf£ (’ Tip path angle = %6.2f degs\n’ ,alphaT*57.3) 
fprinté£ (' Rotor coning angle = %6.2£ degs\n’ ,betao*57.3) 
fprintf(’Location of mean thrust (r/R) = %6.2f \n’,rT2) 
fprintf£ (’ Coliective pitch at .7 r/R = %6.2f degs\n’ , thetao*57.3) 
fprintf (’ solidity = %6.3f \n’,solidity) 
fprintf£ (’ Disk loading = %6.2£ lbs/ft*2\n’ ,DL) 
fprinté (’ Figure of Merit = %6.2£ \n’,FM) 
fprint€ (’ CT/sigma = %6.3f \n’,CT_sig) 
fprintf (' CQ/sigma = %6.4f \n’,CQ sig) 
fprintf (’ CH/sigma = %6.4f \n‘,CH_sig) 
fprintf (' Tip mach of the adv. blade = %6.3f \n’,Machtip) 
fprintf (’ Advance ratio = %6.3f \n’,mu) 
fprintf(’ Rotor thrust required (TPP) = %6.0f lbs\n’,T) 
fprintf (’ Rotor power required = %6.0f h.p.\n’,Protor) 
fprintf (’ Rotor torque = %6.0f ft-1bs\n’ ,Qrotor) 
diary off 
cle 
Gisp(’ ‘') 
disp (’ *** PERFORMANCE OUTPUT DATA INSTRUCTIONS ***’) 
dispi’ ') 
disp (’ A. Single value data saved to a file named:’) 
eval ([‘disp(’’ "' filenamel,’ .prf"’'’)’J) 
disp(’ ’) 
disp (' B. This is a text file, use the TYPE command to view 
the’) 
disp (’ file or use a text editor to view/print the file.’) 
eval (['flagsexist(’’’,filenamel,’ .prf’’);']) 

if flag <1, 
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eval([’!rename diary ’,filenamel,’ .prf’’’)]) 


else 
eval({’!del ',filenamel,’ .prf'’’)) 
eval({’!rename diary ‘',filenamel,’ .prf’’‘)) 
end 


% *** Output to disk (.mat file containing matrix variables 
* xr, psi, vi, theta, betat, alpha, Tpsi, Mpsi, DMpsi, 
% aT, GM, aD) *** 


% *** Configuring variables for output *** 


theta=theta*57.3; 

betat=(betat twist*(0.7- (Reff+(R-Reff) /2) /R)]*57.3; 
alphazalpham*57.3;,alpha=[alpha zeros (psi)); 
Mpsi=Mpsi (:,length(Mpsi(1,:))-1); 

aM= (dM ddM) i 

psi=psi*S7.3; 

r=lr (R-(R-Reff) /2)); 

viz=[vi 0]; 





filename2=[filenamel '_p’]; 


disp(’ ‘) 

disp (’ C. Matrix and vector data saved to a file named:’) 
eval ([’disp(’’ "’ filename2,’ .mat"’’)’)) 

Gisp(’ °) 

disp (’ D. This is a ".mat" binary file, use the LOAD command 
to’) 

disp (’ retrieve the data for plotting.’) 

eval([’save ’,filename2,’ r psi vi theta betat alpha Tpsi Mpsi DMpsi dT 
dM dD’)) 

Gisp(’ ’) 

disp (’ **x* END OF PERFORMANCE ***’) 

disp(’ °) 

disp(’ ‘) 


disp (’ press any key to continue ...’) 








APPENDIX D: THRCALC.M 


sthrealc calculates the total thrust along a blade at 
seach azimuth (psi) location 


Up=zeros (psi*r) ; 
Ut«zeros (Up) ; 
aT=zzeros (Up) ; 
ddTezeros (psi) ; 


for iz1l:length(psi), 


Up (i, :) svi. *cos (betao) +Vinf*sin (alphaT) *cos (betao) +Vinf*cos (alphaT) *sin ( 
betao) *cos (psi (i) ) ; 
Ut (i, :) =x. *omega+Vinf*cos (alphaT) “sin (psi (i) ) ; 
phisatan2 (Up(i,:),Ut(i,:)); 
alphastheta (i) +betat-phi; 
if airfoil==1, 
(CL, CD] shh02clcd (alpha) ; 
else 
{CL, CD] =vr12clced (alpha) ; 
end 


aT (i, :)=0.5*rho*cblade*dr* (Up (i, :) .*24+Ut (i, :) .*2) .*(CL.*cos (phi) -CD.*sin 
(phi) ); 
Tpsi (i) =sum(dT(i,:)); 


%* *** calculations for tip loss area *** 


Uptip=Vinf*sin (alphaT) *cos (betao) +Vinf*cos (alphaT) *sin (betao) *cos (psi (i) 


Uttip=(R- (R-Reff) /2) *omega+Vinf*cos (alphaT) *sin (psi (i) ) ; 
phitipsatan2 (Uptip, Uttip) ; 
ddT (i) =0.5*rho*cblade* (R-Reff) * (Uptip*2+Uttip*2) * (-.009*sin(phitip) ) ; 
Tpsi (i) =Tpsi (i) +ddaT (i) ; 

end 
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APPENDIX E: TMCALC.M 


&tmcalc calculates the total thrust moment along a blade 
wat each azimuth (psi) location 


Up=zeros (psi*r) ; 
Utxzeros (Up) ; 
dM=zeros (Up) ; 
ddM=zeros (psi) ; 


for izl:length (psi), 


Up (i, :) =vi.*cos (betao) +Vinf*sin (alphaT) «cos (betao) +Vinf*cos (alphaT) *sin ( 
betao) *cos (psi (i) ) ; 
Ut (i, :) =r. *omega+Vinf*cos (alphaT) *sin (psi (i)) ; 
phisatan2 (Up(i,:),Ut(i,:)); 
alpha=theta (i,k) +betat-phi; 
if airfoil=<1, 
(CL, CD] =hh02cled (alpha) ; 
else 
(CL, CD] =vr12cled (alpha) ; 
end 


aM (i, :)=0.5*rho*cblade.*r*dr.* (Up (i, :) .*24+Ut (i, :) .*2) .*(CL.*cos (phi) -CD. 
*sin (phi) ) ; 
Mpsi (i,k) =sum(GM(i,:)); 


* *** calculations for tip loss areas *** 


Uptip=Vinf*sin (alphaT) *cos (betao) +Vinf*cos (alphaT) «sin (betao) *cos (psi (i) 


“Uttips (R- (R-Ref£) /2) *omega+Vinf*cos (alphaT) *sin (psi (i)) ; 
phitip=atan2 (Uptip, Uttip) ; 


ddaM (i) =0.5*rho*cblade* (R- (R-Reff) /2) * (R-Reff) * (Uptip*2+Uttip*2) *(-.009*s 
in(phitip)); 

Mpsi (i,k) =Mpsi (i,k) +ddM(i) ; 
end 
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APPENDIX F: DMCALC.M 


tdmcalc calculates the total drag along a blade at 
teach azimuth (psi) location 


Up=zeros (psi*r) ; 
Ut=zeros (Up) ; 
alpham=zeros (Up) ; 
@D=zeros (Up) ; 
ddD=zeros (psi) ; 
ddDM=zeros (psi) ; 


for izl:length (psi), 


Up (i, :) =vi.*cos (betao) +Vinf*sin (alphaT) *cos (betao) +Vinf*cos (alphaT) *sin ( 
betao) *cos (psi (i)); 
Ut (i, :) =r. *omega+Vinf*cos (alphaT) *sin (psi (i) ) ; 
phisatan2 (Up(i,:),Ut(i,:)); 
alphastheta (i) +betat-phi; 
alpham (i, :) =alpha; 
if airfoil== Fi 
(CL, CD} =hhO2clcd (alpha) ; 
else 
(CL, CD) zvr12¢clcd (alpha) ; 
end 


@D (i, :)=0.5*rho*cblade*dr* (Up (i, :) .*2+Ut(i,:).*2).*(CL.*sin (phi) +CD.*cos 
(phi) ) ; 

GDM=dD(i,:) .*r; 

DMpsi (i) =sum (dDM) ; 


* *** calculations for tip loss area *** 


Uptip=Vinf*sin (alphaT) *cos (betao) +Vinf*cos (alphaT) *sin (betao) *cos (psi (i) 


Uttip=(R- (R-Reff) /2) *omega+Vinf*cos (alphaT) *sin (psi (i) ) ; 
phitip=atan2 (Uptip, Uttip) ; 
dadD (i) =0.5*rho*cblade* (R-Reff) * (Uptip*2+Uttip*2) * (.009*cos (phitip) ) ; 
AdGDM (i) =ddD (i) * (R- (R-Reff) /2) ; 
DMpsi (i) =DMpsi (i) +ddDM(i) ; 

end 
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APPENDIX G: HHO2CLCD.M 


function [CL,CD] «nhh02clcd (alpha, Mach) 

thh02clcd calculates CL and CD for an HH-02 airfoil 

tgiven angle of attack (alpha) in radians and Mach number (Mach) 
% [CL, CD) =hh02clcd (alpha, Mach) 

CLezeros (alpha) ; 

CD=zeros (alpha) ; 

asalpha*180/pi; 


chk1=(a><20 & a<=180); 


CL2CL+chki . * (0.42541+0.026863*a+5 .5988e-4*a.°2-2.1493e-S*a.*34+1.5932e-7* 
a.*4-3.4659e-10*a.*S) ; 


CD=CD+chk1 . * (-0.7179+0.061213*a-5.9861e-4*a.~2+7.3708e-6*a.*3-6.6605e-8* 
a.*44+1.913e-10%a.*S) ; 


chki=(a>=-180 & ac=-50); 


CL=CL+chk1.* (-4.6183-0.1923*a-3 .5554e-3*a.*2-3.3273e-5%a.*3-1.4528e-7*a. 
“4-2.3003e-10*a.*5); 


CD=CD+chk1 .* (2.7093e-2-2.1309e-2*a+2.0335e-4*a.~2+3.47e-7%a.*3-3.0586e-8 
*a.*4-1.2584e-10*a.*°5) ; 


chkl=(a>-50 & a<-20); 


CL=CL+chk1.* (-2.5519-0.22847*a-9 .5667e-3*a.°2-1.7051e-4*a.*3-1.0909e-6*a 
“4a); 


CD=CD+chk1 .* (2.7093e-2-2.1309e-2*a+2.0335e-4*a.~2+3.47e-7%*a.*3-3.0586e-8 
wa. *4-1.2584e-10*a.*5) ; 


chkl=(a>=-20 & a<=-10); 
CL=CL+chk1 .* (-0.2+0.089*a+0.0034*a.“*2) ; 


CD=CD+chk1.* (2. 7093e-2-2.1309e-2*a+2.0335e-4*a.~2+3.47e-7*a.~3-3.0586e-8 
*a.*4-1.2584e-10*a.“°5) ; 


chkl=(a<20 & a>-10); 


CL=CL+chk1.* (5.8766e-24+1.313le-1*a+2.4742e-3%a.*2-5.303e-4%a.*3-1.5818e- 
5*a.*44+1.28e-6%a.*5); 

chk2=a<-4; 

ehk2<«chk2.*chki1; 


CD=CD+chk2.* (1.3786+0.916*a+0.21396*a.*2+2.0371le-2*a.*3+7.0076e-4%a.%4) ; 
chk2=(a>=-4 & a<=7); 
chk2=chk2.*chk1; 


CD2CD+chk2 . * (9.732e-3+3 .2326e-4*a+1 .4392e-4%a.*2-8.5073e-5*a.*34+1.1826e- 
6*a.~*441.5271e-6%a.*5); 

chk2<a>7; 

ehk2=chk2.*chk1; 

CD=CD+chk2.* (1.842e-1-5.7532e-2*a+5 .8043e-3*a.*2-1.2803e-4%a.*3); 
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APPENDIX H: VR12CLCD.M 


function [CL,CD] svrl2clcd (alpha) 

‘vrl2cled calculates CL and CD for the VR-12 airfoil 
tgiven angle of attack (alpha) in radians 

% (CL, CD) sevr12cled (alpha) 

CL=zeros (alpha) ; 

CD=zeros (alpha) ; 

azalpha*180/pi; 

chk2z (a>=«20 & a<«180); 


CL«eCL+chk. * (1.1733-0.018879*a+1 .5762e-3%a.*2-3.1925e-5*a.*3+2.0949e-7*a. 
“4-4.3807e-10*a.*5); 


chk= (a>=-180 & a<=-50); 


CLaCL+chk.* (-4.6183-0.1923*a-3.5554e-3*a.*2-3.3273e-S*a.“3-1.4528e-7*a.~ 
4-2.3003e-10*a.“5) ; 


chk=(a>-50 & a<-30); 
CL=CL+chk. * (-0.22114+0.020857*a+2.8571e-4*a.“2); 


chk= (a>=-30 & a<=-10); 
CL2CL+chk.* (-1.11-0.12383*a-0.01515*a.*2-6.8667e-4%a.*3-le-5*a.*4); 


chk=(a<20 & a>-10); 
CL=CL+chk. * (0.11976+0.12341*a+5 .584le-4*a.*2-2.0652e-4*a.*3) ; 


chk=(a>=17 & a<=180) ; 


CD=CD+chk. * (-0.26376+0.017917*a+6 .9927e-4%a. ~*2-9.1137e-6%a.~3+2.6277e-8* 
a.*4); 


chk=(a>=-180 & ac=-10); 


CD=CD+chk.* (-0.17486-0.034463*a-1.0233e-4%*a.*2-2.8958e-6%a.“3-4.6577e-8* 
a.*4-1.5557e-10*a.°S); 


chk=(a>-10 & a<=0); 


CD=CD+chk. * (9.8678e-3+3 .4934e-3*a+1.4844e-3*a.*2-1.3564e-4%a.*3-1.0936e- 
5*a.*4); 


chk=(a>0 & a<=15); 


CD=CD+chk.* (9. 8e-3+7.0457e-4*a+5 .6104e-5*a.*2-4.1151le-5*a.*3+3.8695e-6%a 
.*4); 


chk=(a>15 & a<i7); 
CD=CD+chk.* (-1.33+1.325e-1*a-2.5e-3*a.*2) ; 
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